---
title: "Effects of obstruent voicing on vowel fundamental frequency in Dutch: Supplementary Materials"
author: "Hugo Quené & Anne-France Pinget"
date: "8 March 2023"
output: pdf_document
---

These supplementary materials provide further details and additional outputs of the GAMMs reported in our manuscript for _J. Acoust. Soc. Am._ 

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
setwd("~/SURFdrive/onderzoek/annefrancepinget/SupplementaryMaterial")
options(repos="https://CRAN.R-project.org") # override user setting
```

# Load and clean data

## Load prepared data

We start the analysis by loading the pre-processed data set `carrier2B194v2`. This data set was prepared according to script `AFPHQprepare20221123a.Rmd`.

We also set `P50`, the median duration (in ms) of vowel tokens, for later use in figures.  

```{r save.or.load.data}
load( file="carrier2B194v2.Rdata" ) # object name is carrier2B
P50 <- 97 # ms, lower weights after this point
```

> This dataset has 49342 observations, from 4492 vowel tokens, with 11 observations from each vowel token.  

## Ignore suspect f0 jumps

As dependent variable we use `measure.c4`, which is the f0 centered to the speaker's centroid, in semitones, and with that centroid based on the first 20 ms of the vowel token. 

There are a few observations with very large f0 jumps, defined here as consecutive measurements within the same vowel token with an absolute difference of >6 semitones relative to the previous measurement in the same vowel token. 
These are likely to have been due to f0 measurement errors, and not to true f0 contours. Because these have resulted in unstable GAMMS (i.e. to overfitting of the smoothed f0 contours in the GAMM), we will ignore these few suspect measurements. 

Column `start_event` is a boolean or logical variable, indicating whether this observation is the first f0 observation of a new vowel token (event).

```{r f0jumps}
with( carrier2B194v2, diff(measure.c4) ) -> diff.measure.c4
diff.measure.c4 <- c(NA,diff.measure.c4)
which( abs(diff.measure.c4)>6 & !carrier2B194v2$start.event) -> jumps
# manual inspection here: suspect row numbers in carrier2B194v2
f0jumps <- c( 7235:7236,
              15081,
              18239:18240,
              40588:40598,
              41192,
              47868,
              48362:48363 )
length(f0jumps)
( unique(carrier2B194v2$file[f0jumps]) )
# rm(diff.measure.c4,jumps)
```

This results in ignoring 20 f0 observations from 7 vowel tokens. 

## Set color palette 

For colour schemes used below, see 
https://personal.sron.nl/~pault/data/colourschemes.pdf.
Here we use a *vibrant* color scheme (Fig.3 in the linked pdf), with warm colors for fricatives and cool colors for plosives. Lighter for voiced and darker for unvoiced consonants.  
Colours are specified directly, so as to avoid installing additional packages containing palette definitions, etc.  

```{r set-palette}
# make four vibrant colours
cp <- rgb( red  =c(000,000,238,204)/255,
           green=c(119,153,119,051)/255,
           blue =c(187,136,051,017)/255,
           names=c("blue","teal","orange","red") )
cp <- c( "black", cp[c(2,4,1,3)]) # as by consonant levels: m b f p v
```

---

# Model 1: general GAMM over all consonants

```{r gamm-setup}
require(mgcv)
require(itsadug)
```

## Define GAMM

First we define the ideal full model. This model includes `cons` (the identity of the preceding consonant), smoothed effect of `vowtime`, and the interactions. All terms are allowed to vary over participants (random slopes). 

Here we use a different smoothed curve for each consonant, using the `by` argument in the `s` call; see https://www.r-bloggers.com/2017/07/generalized-addictive-models-and-mixed-effects-in-agriculture/ and see Wieling 2018 (J. Phon.). 

As suggested by pilot analyses, we use a `scat` distribution (not gaussian but scaled-t distribution) of residuals, with wider tails.


## Run GAMM 1

Using suggestions and ideas from Wieling (2018, J Phon), in particular his models `m4` and m5` in sections 4.7.3 and 4.7.3.

The data set with suffix `v2` is sorted by file and vowtime, and it contains a column `start_event`. 

The GAMM is already estimated, and during knitting this R notebook, the GAMM is read from disk. 

```{r m1s4-gamm}
set.seed(20230307) # for replicability
# commented out for knitting
# system.time(
#   m1s4.gamm <- bam( measure.c4 ~ s(vowtime,by=cons,k=20) + cons +
#                        ti(vowtime,ppn,bs="fs",m=1), # re slopes
#                   rho=.85, # from pilot models
#                   gamma=1.6, # Wood, 2017, p.186
#                   AR.start=start.event,
#                   family=scat(min.df=3), # scaled t distrib
#                   # nthreads = 2, # 2 cores
#                   discrete=TRUE,
#                   data=carrier2B194v2[-f0jumps,], # ignoring f0 jumps
#                   weights=wght ) # using weights
# ) # ~42 seconds, no warnings
# save( m1s4.gamm, file="./m1s4gamm20230307.Rdata" ) # for knitting
load( file="./m1s4gamm20230307.Rdata" ) # for knitting
```

Check the first GAMM...

```{r m1-check-one}
gam.check(m1s4.gamm) # no problems with edf; k-index<1 but p looks good
summary( m1s4.gamm )
m1s4.gamm.acf <- acf_resid(m1s4.gamm) # |rho|<.1, good
# gamtabs(m1s4.gamm)
```

Rank check is OK. 

The smoothing functions require about 5 to 15 knots. (Although this may suggest possible overfitting of the f0 contour, over 11 points per vowel token, models with lower `k` did not pass the check of `edf` being substantially lower than $k-1$).

Deviance explained is 38.4 % (which is higher than with `family=gaussian` GAMM during piloting).  

The residuals look good, and autocorrelation is very low. 

The summary data are reported in Appendix Table I of the manuscript. 

## Figure 1

```{r m1s-check-two}
# load(file="./m1s4.gamm.Rdata")
require(itsadug)
# P50 <- 97 # added 2023.01.11 for completeness; P50 of vowel durations
plot_smooth( m1s4.gamm, 
             view="vowtime", plot_all="cons", 
             rug=FALSE, col=cp, lwd=3, 
             lty=c(1,5,4,3,6), # m b f p v
             legend_plot_all = c(x=50,y=0),
             xlim=c(0,P50*1.04), # allow plot room
             xlab="vowel time (ms)", ylab="f0 (ST, centered)")
legend( x=c(85,100), y=c(0,-1.5), col=cp, lwd=3, bty="n", 
        lty=c(1,5,4,3,6), # m b f p v
        legend=levels(carrier2B194v2$cons) ) -> fig1legend
abline(v=50,col="grey50",lty=3)
pdf( file="./m1s4gamm.v20230307.pdf", width=6, height=4)
plot_smooth( m1s4.gamm,
             view="vowtime", plot_all="cons",
             rug=FALSE, col=cp, lwd=3,
             lty=c(1,5,4,3,6), # m b f p v
             xlim=c(0,P50*1.04), # allow plot room
             xlab="vowel time (ms)", ylab="f0 (ST, centered)")
legend( x=c(85,100), y=c(0,-1.5), col=cp, lwd=3, bty="n", 
        lty=c(1,5,4,3,6), # m b f p v
        legend=levels(carrier2B194v2$cons) ) -> fig1legend
dev.off()
```
Fig. 1. (Color online). Fitted contours of f0 (in semitones, centered by speaker) over time (in ms from vowel onset) as fitted by nonlinear smooth functions, broken down by preceding consonant, with 95% confidence intervals. Random effects are ignored in the fitted f0 values. [file m1s4gamm.v20230307.pdf].

> Figure 1 shows several interesting patterns. Firstly, compared to the baseline condition of f0 after /m/, the f0 after unvoiced consonants /f, p/ starts relatively higher at vowel onset (by about 2.0 and 1.5 semitones, respectively, see Figure 1). In contrast, secondly, the f0 after voiced consonants /v, b/ does not differ from that after /m/. At vowel onset, the confidence bands of f0 after unvoiced consonants /b, p/ are at least 1 semitone higher than those after voiced /v, f/. Third, these voicing effects of the preceding consonant on f0 are only transitory, and are no longer noticeable after approximately 50 ms in the vowel (when the 95% confidence intervals of /m/ and /p/ start to overlap, see Figure 1). 

---

# Model 2: after fricatives

For fricatives, we include the degree of voicing of the preceding consonant as a predictor. This requires some data wrangling...

## Wrangle voicing data

The predictor `cons` denoting the preceding word-initial consonant is a factor with `m` as baseline.  

```{r levels-of-variable}
levels(carrier2B194v2$cons) # triple check
```

Next, we include the degree of voicing of the preceding consonant (`voicing_A`) as a variable in the data set. We do this for _all_ preceding consonants. For /m/ the raw degree of voicing is always 100%. 

The degree of voicing itself is inspected first, by taking only a single observation for each vowel token, at vowel _midpoint_ (where `type_measure==0.5`). (This turns out to be safer than using the vowel _onset_ where `type_measure==0`).

```{r consonants}
carrier2consA <- subset( carrier2B194v2, # for cons data
                         type_measure==0.5 )
```

## Summary of raw `voicing_A` 

Each of the five preceding consonants in this study occurred in 9 tokens (see Appendix of manuscript) spoken by 100 speakers. 

```{r cons-voicing-summary-2}
require(dplyr)
carrier2consA %>% 
  group_by(cons) %>% 
  summarize(mean=mean(voicing_A), sd=sd(voicing_A), n=length(voicing_A) ) 
carrier2plosA <- subset( carrier2consA, cons %in% c("p","b") )
carrier2fricA <- subset( carrier2consA, cons %in% c("f","v") )
```

Note that for the fricatives' raw degree of voicing, SD for /v/ is considerably higher than for /f/. 

## Centering of `voicing_A`

First we center the raw `voicing_A` predictor. 

Here we use the max voicing value, 100%, as the reference value for centering. This means that the max voicing value of 100%, as found for /m/, is centered to zero. The polarity of the scaled voicing is preserved: higher values denote more voicing. (This centering is an artefact of other unreported analyses.)

> For easier interpretation, the centered scale from $-100$ to 0 will be uncentered to run from 0 to 100 in Figure 2. 

Because the `voicing_A` parameter is bimodally distributed, no scaling is done. 
The centered values are stored as a separate column `voi.c` in a new data frame, with suffix `v3`, which is also saved.

```{r cons-voicing-scale}
carrier2B194v3 <- cbind( carrier2B194v2,
                         voi.c = scale(carrier2B194v2$voicing_A,
                               center=100, scale=FALSE) )
save( carrier2B194v3, file="./carrier2B194v3.Rdata" )
```

We extract two data sets for the two consonant categories, with a single observation per vowel token (at `type_measure`=0.5, i.e., at vowel midpoint). These datasets are named `A`, and these contain consonant data (a single row per token). 

We also extract two similar data sets, for the two consonant categories, with all observations per vowel token. These datasets are named `B`, and these contain vowel data (multiple rows per token).

```{r cons-voicing-scale-A}
carrier2plosA <- subset( carrier2B194v3, 
                         (cons %in% c("p","b")) & (type_measure==0.5) )
carrier2plosB <- subset( carrier2B194v3, 
                         (cons %in% c("p","b")) )
carrier2fricA <- subset( carrier2B194v3, 
                         (cons %in% c("f","v")) & (type_measure==0.5) )
carrier2fricB <- subset( carrier2B194v3, 
                         (cons %in% c("f","v")) )
```

For consistency, we re-use a previously prepared dataset for vowels after fricatives... 

```{r fric-BB-dataset}
load( file="./carrier2fricBB.Rdata" ) # fric+V f0 data
```

## Run GAMM 2

... and we use the same random seed as used in the analysis reported in the manuscript. Here we load the previously estimated second GAMM from disk, rather than computing it again (which would take >5 hours). 

```{r gamm2}
## voicing_A, raw, pooled over consonants, carrier2fricBB
set.seed(20221127) # for replicability
# commented out for knitting
# print( paste("Estimating fricative models m100 and m101 started at",Sys.time()) )
# system.time( # null
#   fric.m100fs3.gamm <- bam( measure.c4 ~ s(vowtime,by=cons) + cons +
#                              s(vowtime,ppn,by=cons,bs="fs",m=2) + # random slopes
#                              s(voi.c,by=cons,k=20), # voicing, raw
#                            AR.start=start.event,
#                            discrete=TRUE, # for speeding up
#                            family="scat",
#                            # nthreads = 4, # 4 cores
#                            data=carrier2fricB, weights=wght ) # BB not A
# ) # 624 seconds; null model, warnings ignored
# system.time( # full
#   fric.m101fs3.gamm <- bam( measure.c4 ~ cons +
#                              te(voi.c,vowtime,by=cons,k=20) +
#                              s(vowtime,ppn,by=cons,bs="fs",m=2), # random slopes
#                            AR.start=start.event,
#                            discrete=TRUE, # for speeding up
#                            family="scat",
#                            # nthreads = 4, # 4 cores
#                            data=carrier2fricBB, weights=wght ) # B not A
# ) # 19201 seconds; full model, no warnings
# print( paste("Finished at",Sys.time()) )
# save( fric.m100fs3.gamm,
#       file="./fric.m100fs3.gamm.Rdata" ) # comment for knitting
# save( fric.m101fs3.gamm,
#       file="./fric.m101fs3.gamm.Rdata" ) # comment for knitting
load( file="./fric.m100fs3.gamm.Rdata" ) # uncomment for knitting
load( file="./fric.m101fs3.gamm.Rdata" ) # uncomment for knitting
```

Check the second GAMM...

```{r m2-check}
gam.check(fric.m101fs3.gamm) 
print( summary( fric.m101fs3.gamm ) -> fric.m101fs3.gamm.summary )
anova( fric.m100fs3.gamm, fric.m101fs3.gamm, test="Chisq")
```


Rank check is OK. 

Two of the k-indices have $p<.05$, but the edf are always well below $k-1$. 

Deviance explained is 51.1 % .

The model containing the interaction effect provides a far better fit than the model without, $\chi^2(73)=385, p<.001$.

## Figure 2

This is a two-panel figure, left after /v/, right after /f/.
We also need the values of `voi.c` from the fricative consonant dataset `carrier2fricAA`. 

We suppress value labels along the centered (shifted) axis of degree of voicing, and add axis labels later. 

```{r figure2}
par( mfcol=c(1,2)) -> op # two panels
require(vioplot) # add half violin plot along Y axis
# /v/ panel
thecons <- "v"
load(file="./carrier2fricAA.Rdata") # consonant data, already loaded when knitting
ymarks <- carrier2fricAA[carrier2fricAA$cons==thecons, "voi.c"]
fvisgam( fric.m101fs3.gamm,
         view=c("vowtime","voi.c"),
         cond=list(cons=c(thecons)),
         xlab="vowel time (ms)",
         yaxt="n", 
         # ylab=paste("pct. voicing of /",thecons,"/",sep=""),
         ylab="",
         contour.col="black", # not cp[5] # v is level 5 of cons
         color="gray",
         xlim=c(2,P50*.8), # P50 of vowel duration
         ylim=c(-97,-3), # range of voi.c of /v/
          zlim=c(-1,2),
         show.diff=TRUE, col.diff="white", alpha.diff=.99,
         main=paste("after /",thecons,"/",sep=""),
         add.color.legend = FALSE,
         n.grid=15,
         rm.ranef=TRUE ) -> fricv.fvisgam
# abline( h=cons.pc1.medians[thecons], lty="dashed", lwd=2, col=cp[3] )
axis( side=2, at=seq(-100,+0,by=10), 
      labels=seq(-100,0,by=10)+100 )
mtext( paste("percentage voicing of /",thecons,"/", sep=""), side=2, line=2.5 )
# rug( ymarks, side=2, ticksize=0.04, quiet=TRUE, col="black" )
vioplot( I(ymarks), at=-0.5, xaxt="n", side="right", add=T, 
         # yaxt="n", ylim=c(-95,-5), # range of voi.c of /v/
         h=h.select(ymarks,method="cv"), # h=6 
         col="grey50",
         plotCentre="point", wex=14 )
# /f/ panel
thecons <- "f"
ymarks <- carrier2fricA[carrier2fricA$cons==thecons, "voi.c"]
fvisgam( fric.m101fs3.gamm,
         view=c("vowtime","voi.c"),
         cond=list(cons=c(thecons)),
         xlab="vowel time (ms)",
         yaxt="n", 
         # ylab=paste("pct. voicing of /",thecons,"/",sep=""),
         ylab="",
         contour.col="black", # not cp[5] # v is level 5 of cons
         color="gray",
         xlim=c(2,P50*.8), # P50 of vowel duration
         ylim=c(-98,-65), # range of voi.c of /v/
         zlim=c(-1,2),
         show.diff=TRUE, col.diff="white", alpha.diff=.99,
         main=paste("after /",thecons,"/",sep=""),
         add.color.legend = FALSE,
         n.grid=15,
         rm.ranef=TRUE ) -> fricf.fvisgam
# abline( h=cons.pc1.medians[thecons], lty="dashed", lwd=2, col=cp[3] )
axis( side=2, at=seq(-90,0,by=10), 
      labels=seq(-90,0,by=10)+100 )
mtext( paste("percentage voicing of /",thecons,"/", sep=""), side=2, line=2.5 )
# rug( ymarks, side=2, ticksize=0.04, quiet=TRUE, col="black" )
# h.select(ymarks,method="cv") # CV yields h=6 as smoothing param
vioplot( I(ymarks), at=-0.5, xaxt="n", side="right", add=T, 
         h=h.select(ymarks,method="cv"), # h=6 
         col="grey50",
         plotCentre="point", wex=14*1.5 )
par(op) # restore earlier graphics settings
```

Fig. 2. Contour plots showing the interaction between vowel time (in ms, first part of vowel only) and the degree of voicing of the preceding consonant (in percentage points, rescaled), on f0 (in semitones, centered, lighter shade indicates higher f0). The lefthand panel shows the interaction for vowel tokens after /v/; the righthand panel shows the interaction after /f/; voicing axes are scaled differently across panels for clarity. The density of the degree of voicing is shown along the left axis (with quartiles and median). Random effects have been ignored in the fitted f0 values. White areas have an essentially flat contour. [file fric101gamm20230303.pdf]

## Supplementary figure

This figure highlights the difference in distributions of the degrees of voicing of /v/ and /f/.  

```{r voi.distributions}
require(vioplot)
voidata <- subset(carrier2fricA, cons %in% c("v","f") )
voidata$cons <- as.character(voidata$cons) # to exclude missing levels of cons
vioplot( voi.c~cons, data=voidata,
         xlab="Preceding consonant",
         ylab="Degree of voicing (centered to 100)" )
```

---

# Model 3: after plosives

For consistency (and as with fricative consonants), we re-use a previously prepared dataset for vowels after plosives... 

```{r plos-BB-dataset}
load( file="./carrier2plosBB.Rdata" ) # plos+V f0 data
```

## Run GAMM 3

... and we use the same random seed as used in the analysis reported in the manuscript. Here we load the previously estimated second GAMM from disk, rather than computing it again (which would take >5 hours). 

```{r gamm3}
## VOT_ms, raw, pooled over consonants, carrier2plosBB
set.seed(20221127) # for replicability
# commented out for knitting
# print( paste("Estimating plosive models m100 and m101 started at",Sys.time()) )
# system.time( # null
#   plos.m100fs3.gamm <- bam( measure.c4 ~ s(vowtime,by=cons) + cons +
#                              s(vowtime,ppn,by=cons,bs="fs",m=2) + # random slopes
#                              s(VOT_ms,by=cons,k=20), # VOT, raw
#                            AR.start=start.event,
#                            discrete=TRUE, # for speeding up
#                            family="scat", nthreads = 4, # 4 cores
#                            data=carrier2plosBB, weights=wght ) # BB not AA
# ) # 
# system.time( # full
#   plos.m101fs3.gamm <- bam( measure.c4 ~ cons +
#                              te(VOT_ms,vowtime,by=cons,k=20) +
#                              s(vowtime,ppn,by=cons,bs="fs",m=2), # random slopes
#                            AR.start=start.event,
#                            discrete=TRUE, # for speeding up
#                            family="scat", nthreads = 4, # 4 cores
#                            data=carrier2plosBB, weights=wght ) # BB not AA
# ) # seconds
# print( paste("Finished at",Sys.time()) )
# save( plos.m100fs3.gamm, 
#       file="./plos.m100fs3.gamm.Rdata" ) # comment for knitting
# save( plos.m101fs3.gamm, 
#       file="./plos.m101fs3.gamm.Rdata" ) # comment for knitting
# do not load if models are still in global environment
set.seed(20221129) # for replicability
load( file="./plos.m100fs3.gamm.Rdata" ) # comment for knitting
load( file="./plos.m101fs3.gamm.Rdata" ) # comment for knitting
```

Check model ...

```{r m3-check}
gam.check(plos.m101fs3.gamm) # good
print( summary( plos.m101fs3.gamm ) -> plos.m101fs3.gamm.summary )
anova( plos.m100fs3.gamm, plos.m101fs3.gamm, test="Chisq") # p=.06
```


Rank check is OK. 

The k-indices are also OK (>1). 

Deviance explained is 49.9 % .

The model containing the interaction effect does not provide a significantly better fit than the model without, $\chi^2(48)=63.3, p=.063$.

## Figure 3

This is a two-panel figure, left after /b/, right after /p/.
We also need the values of raw `VOT_ms` (in ms) from the plosive consonant dataset `carrier2plosAA`. 

```{r figure3}
par( mfcol=c(1,2)) -> op # two panels
require(vioplot) # add half violin plot along Y axis
# /b/ panel
thecons <- "b"
load("./carrier2plosAA.Rdata") # consonant data, already loaded when knitting
ymarks <- carrier2plosAA[carrier2plosAA$cons==thecons, "VOT_ms"]
fvisgam( plos.m101fs3.gamm,
         view=c("vowtime","VOT_ms"),
         cond=list(cons=c(thecons)),
         xlab="vowel time (ms)",
         ylab=paste("VOT of /",thecons,"/ (ms)",sep=""),
         xlim=c(3,P50*.8), # P50 of vowel duration
         ylim=c(-150,-50), # range of VOT of plosives
         zlim=c(-3,3), 
         show.diff=TRUE, col.diff="white", alpha.diff=.99,
         main=paste("after /",thecons,"/",sep=""),
         add.color.legend = FALSE,
         n.grid=12,
         rm.ranef=TRUE, color="gray" ) -> plosb.fvisgam
vioplot( I(ymarks), at=-0.25, xaxt="n", side="right", add=T, 
         h=h.select(ymarks,method="cv"), # h=6 
         col="grey50",
         plotCentre="point", wex=14*1.5 )
# /p/ panel
thecons <- "p"
ymarks <- carrier2plosA[carrier2plosA$cons==thecons, "VOT_ms"]
fvisgam( plos.m101fs3.gamm,
         view=c("vowtime","VOT_ms"),
         cond=list(cons=c(thecons)),
         xlab="vowel time (ms)",
         ylab=paste("VOT of /",thecons,"/ (ms)",sep=""),
         #          contour.col=cp[5], # v is level 5 of cons
         #          color="gray",
         xlim=c(3,P50*.8), # using P50 of vowel duration
         ylim=c(5,45), # range of VOT of plosives
         zlim=c(-2.5,2.5), # range of semitone scale
         show.diff=TRUE, col.diff="white", alpha.diff=.99,
         main=paste("after /",thecons,"/",sep=""),
         add.color.legend = FALSE,
         n.grid=12,
         rm.ranef=TRUE, color="gray" ) -> plosp.fvisgam
vioplot( I(ymarks), at=-0.25, xaxt="n", side="right", add=T, 
         h=h.select(ymarks,method="cv"), # h=6 
         col="grey50",
         plotCentre="point", wex=14*1.5 )
par(op) # restore earlier graphics settings
```

Fig. 3. Contour plots showing the interaction between vowel time (in ms, first part of vowel only) and voice onset time (VOT, in ms), on f0 (in semitones, centered, lighter shade indicates higher f0). The lefthand panel shows the interaction for vowel tokens after /b/; the righthand panel shows the interaction after /p/; VOT axes are scaled differently across panels for clarity. The density of VOTs is shown along the left axis (with quartiles and median). Random effects have been ignored in the fitted f0 values. White areas have an essentially flat contour. [file plos101gamm20230303.pdf]

---

# Cleanup

```{r cleanup}
rm( fric.m100fs3.gamm, plos.m100fs3.gamm )
rm( op )
```

---

# Admin

```{r admin}
sessionInfo()
citation() # base R
citation("mgcv")
citation("itsadug")
```

---

HQ <h.quene@uu.nl> 2023-03-08
