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

Packages needed:

library(ggplot2)
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(car)
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
library(xtable)
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.6.3
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## Warning: package 'tidyr' was built under R version 3.6.2
## 
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack
library(texreg)
## Version:  1.36.23
## Date:     2017-03-03
## Author:   Philip Leifeld (University of Glasgow)
## 
## Please cite the JSS article in your publications -- see citation("texreg").
## 
## Attaching package: 'texreg'
## The following object is masked from 'package:tidyr':
## 
##     extract
library(MuMIn)
## Registered S3 method overwritten by 'MuMIn':
##   method         from
##   predict.merMod lme4
library(ggpubr)
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:texreg':
## 
##     extract
## The following object is masked from 'package:tidyr':
## 
##     extract
library(knitr)
## Warning: package 'knitr' was built under R version 3.6.2
library(emmeans)
## Warning: package 'emmeans' was built under R version 3.6.1

1 Acoustic analysis

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


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

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

1.1 Data used

Unequal data points for respective target words are due to annotation decisions

Monosyllabic

table(mono[mono$position=="phrase-medial",]$word,mono[mono$position=="phrase-medial",]$sylpos)
##        
##          C  N  O
##   Bahn  47 47 47
##   Bann  46 46 46
##   Beet  47 47 47
##   Bett  49 49 49
##   Ruhm  50 50 50
##   Rum   50 50 50
##   Stiel 48 48 48
##   still 50 50 50
table(mono[mono$position=="phrase-final",]$word,mono[mono$position=="phrase-final",]$sylpos)
##        
##          C  N  O
##   Bahn  50 50 50
##   Bann  52 52 52
##   Beet  51 51 51
##   Bett  50 50 50
##   Ruhm  48 48 48
##   Rum   49 49 49
##   Stiel 49 49 49
##   still 48 48 48
addmargins(table(mono$position,mono$sylpos))
##                
##                    C    N    O  Sum
##   phrase-medial  387  387  387 1161
##   phrase-final   397  397  397 1191
##   Sum            784  784  784 2352
addmargins(table(mono$word,mono$sylpos))
##        
##            C    N    O  Sum
##   Bahn    97   97   97  291
##   Bann    98   98   98  294
##   Beet    98   98   98  294
##   Bett    99   99   99  297
##   Ruhm    98   98   98  294
##   Rum     99   99   99  297
##   Stiel   97   97   97  291
##   still   98   98   98  294
##   Sum    784  784  784 2352
m<-addmargins(table(mono$word,mono$sylpos))

784 target words were analyzed in total

Disyllabic

table(di[di$position=="phrase-medial",]$word,di[di$position=="phrase-medial",]$sylpos)
##         
##          N1 N2 O1 O2
##   Huete  48 48 48 48
##   Huette 50 50 50 50
##   Miete  49 48 49 49
##   Mitte  48 48 48 48
table(di[di$position=="phrase-final",]$word,di[di$position=="phrase-final",]$sylpos)
##         
##          N1 N2 O1 O2
##   Huete  50 50 50 50
##   Huette 48 48 48 48
##   Miete  49 49 49 49
##   Mitte  50 50 50 50

197 target words were analyzed

addmargins(table(di$position,di$sylpos))
##                
##                   N1   N2   O1   O2  Sum
##   phrase-medial  195  194  195  195  779
##   phrase-final   197  197  197  197  788
##   Sum            392  391  392  392 1567
addmargins(table(di$word,di$sylpos))
##         
##            N1   N2   O1   O2  Sum
##   Huete    98   98   98   98  392
##   Huette   98   98   98   98  392
##   Miete    98   97   98   98  391
##   Mitte    98   98   98   98  392
##   Sum     392  391  392  392 1567
addmargins(table(di$word,di$position))
##         
##          phrase-medial phrase-final  Sum
##   Huete            192          200  392
##   Huette           200          192  392
##   Miete            195          196  391
##   Mitte            192          200  392
##   Sum              779          788 1567
d<-addmargins(table(di$word,di$sylpos))

Total data excluded

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

1.2 Analysis of monosyllabic target words

### Analysis ####

Normalisation to speech rate: sylposdur (ms) divided by 1000 (sec)., multiplied with (sil/s), results in syllables

meta<-read.table("meta.speech.dur.txt",sep="\t", fileEncoding = "UTF-8")
meta$subject
##  [1] f1 f2 f3 m1 f4 m2 m3 f5 m4 m5
## Levels: f1 f2 f3 f4 f5 m1 m2 m3 m4 m5
# merge data with speech rate
mono<-merge(mono, meta[,c("subject","sprechgeschw")], by = "subject")
di<-merge(di, meta[,c("subject","sprechgeschw")], by = "subject")

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

head(as.data.frame(mono))
##   subject labels    start      end                                        utts
## 1      f1      b 2210.802 2238.177 D_vp01_w1_wav_tg_output_bdl:DP1_Bahn_1_0005
## 2      f1     a: 2243.177 2376.802 D_vp01_w1_wav_tg_output_bdl:DP1_Bahn_1_0005
## 3      f1      n 2376.802 2433.385 D_vp01_w1_wav_tg_output_bdl:DP1_Bahn_1_0005
## 4      f1      b 1796.260 1832.594 D_vp01_w1_wav_tg_output_bdl:DP1_Bahn_2_0074
## 5      f1     a: 1839.573 1967.844 D_vp01_w1_wav_tg_output_bdl:DP1_Bahn_2_0074
## 6      f1      n 1967.844 2010.781 D_vp01_w1_wav_tg_output_bdl:DP1_Bahn_2_0074
##                                db_uuid                     session
## 1 1dbcb4ee-dcc2-11e6-9ce0-e15c0fc30b38 D_vp01_w1_wav_tg_output_bdl
## 2 1dbcb4ee-dcc2-11e6-9ce0-e15c0fc30b38 D_vp01_w1_wav_tg_output_bdl
## 3 1dbcb4ee-dcc2-11e6-9ce0-e15c0fc30b38 D_vp01_w1_wav_tg_output_bdl
## 4 1dbcb4ee-dcc2-11e6-9ce0-e15c0fc30b38 D_vp01_w1_wav_tg_output_bdl
## 5 1dbcb4ee-dcc2-11e6-9ce0-e15c0fc30b38 D_vp01_w1_wav_tg_output_bdl
## 6 1dbcb4ee-dcc2-11e6-9ce0-e15c0fc30b38 D_vp01_w1_wav_tg_output_bdl
##            bundle start_item_id end_item_id level start_item_seq_idx
## 1 DP1_Bahn_1_0005            43          43   MAU                 13
## 2 DP1_Bahn_1_0005            45          45   MAU                 15
## 3 DP1_Bahn_1_0005            46          46   MAU                 16
## 4 DP1_Bahn_2_0074            43          43   MAU                 13
## 5 DP1_Bahn_2_0074            45          45   MAU                 15
## 6 DP1_Bahn_2_0074            46          46   MAU                 16
##   end_item_seq_idx    type sample_start sample_end sample_rate      position
## 1               13 SEGMENT       106119     107432       48000 phrase-medial
## 2               15 SEGMENT       107673     114086       48000 phrase-medial
## 3               16 SEGMENT       114087     116802       48000 phrase-medial
## 4               13 SEGMENT        86221      87964       48000 phrase-medial
## 5               15 SEGMENT        88300      94456       48000 phrase-medial
## 6               16 SEGMENT        94457      96517       48000 phrase-medial
##   word iteration pausedur labelsdur labelsfac sylpos tenseness sylposdur
## 1 Bahn         1        0  27.37500         b      O     tense  32.37500
## 2 Bahn         1        0 133.62500        a:      N     tense 133.62500
## 3 Bahn         1        0  56.58333         n      C     tense  56.58333
## 4 Bahn         2        0  36.33333         b      O     tense  43.31250
## 5 Bahn         2        0 128.27083        a:      N     tense 128.27083
## 6 Bahn         2        0  42.93750         n      C     tense  42.93750
##   sylseglabel sprechgeschw sylposdurnorm
## 1          bh         5.53     0.1790338
## 2          a:         5.53     0.7389463
## 3           n         5.53     0.3129058
## 4          bh         5.53     0.2395181
## 5          a:         5.53     0.7093377
## 6           n         5.53     0.2374444
mono$sylpos<-factor(mono$sylpos, levels=c("O","N","C"))

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

# Overall model monosyllabic acoustic ####
m1<-lmer(sylposdur ~ position + sylpos + tenseness + (1|subject) + (1|word), data = mono)
plot(allEffects(m1))

qqPlot(resid(m1))

## [1] 1890 1891
m2<-lmer(sylposdur ~ position * sylpos + tenseness + (1|subject) + (1|word), data = mono)
anova(m1,m2)
## refitting model(s) with ML (instead of REML)
## Data: mono
## Models:
## m1: sylposdur ~ position + sylpos + tenseness + (1 | subject) + (1 | 
## m1:     word)
## m2: sylposdur ~ position * sylpos + tenseness + (1 | subject) + (1 | 
## m2:     word)
##    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m1  8 23784 23830 -11884    23768                             
## m2 10 23529 23587 -11755    23509 258.82      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word), data = mono)
anova(m2,m3)
## refitting model(s) with ML (instead of REML)
## Data: mono
## Models:
## m2: sylposdur ~ position * sylpos + tenseness + (1 | subject) + (1 | 
## m2:     word)
## m3: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 | 
## m3:     word)
##    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m2 10 23529 23587 -11755    23509                             
## m3 15 23222 23308 -11596    23192 317.27      5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m4<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+position|subject) + (0+position|word), data = mono)
## boundary (singular) fit: see ?isSingular
anova(m3,m4,refit=F)
## Data: mono
## Models:
## m3: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 | 
## m3:     word)
## m4: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 | 
## m4:     word) + (0 + position | subject) + (0 + position | word)
##    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m3 15 23172 23258 -11571    23142                             
## m4 21 23152 23273 -11555    23110 31.364      6   2.16e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(m4))

summary(m4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 |  
##     word) + (0 + position | subject) + (0 + position | word)
##    Data: mono
## 
## REML criterion at convergence: 23110.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7961 -0.6901 -0.0360  0.5685  4.3391 
## 
## Random effects:
##  Groups    Name                  Variance  Std.Dev.  Corr 
##  subject   positionphrase-medial 5.288e+01  7.272055      
##            positionphrase-final  2.326e+01  4.823273 -0.13
##  subject.1 (Intercept)           5.164e+01  7.186317      
##  word      positionphrase-medial 1.305e+02 11.424254      
##            positionphrase-final  2.306e+02 15.185981 1.00 
##  word.1    (Intercept)           3.890e-06  0.001972      
##  Residual                        1.079e+03 32.850649      
## Number of obs: 2352, groups:  subject, 10; word, 8
## 
## Fixed effects:
##                                           Estimate Std. Error       df t value
## (Intercept)                                 96.426      6.979   11.205  13.817
## positionphrase-final                         9.897      4.808   29.080   2.059
## sylposN                                      3.725      3.353 2315.958   1.111
## sylposC                                    -31.967      3.353 2315.958  -9.534
## tensenesslax                                -7.907      8.742    7.396  -0.905
## positionphrase-final:sylposN                34.814      4.706 2315.958   7.399
## positionphrase-final:sylposC                53.826      4.706 2315.958  11.439
## positionphrase-final:tensenesslax           -2.155      5.396   32.527  -0.399
## sylposN:tensenesslax                       -29.446      4.723 2315.958  -6.234
## sylposC:tensenesslax                        15.801      4.723 2315.958   3.345
## positionphrase-final:sylposN:tensenesslax  -16.700      6.637 2315.958  -2.516
## positionphrase-final:sylposC:tensenesslax   10.663      6.637 2315.958   1.606
##                                           Pr(>|t|)    
## (Intercept)                               2.19e-08 ***
## positionphrase-final                      0.048595 *  
## sylposN                                   0.266632    
## sylposC                                    < 2e-16 ***
## tensenesslax                              0.394228    
## positionphrase-final:sylposN              1.92e-13 ***
## positionphrase-final:sylposC               < 2e-16 ***
## positionphrase-final:tensenesslax         0.692178    
## sylposN:tensenesslax                      5.38e-10 ***
## sylposC:tensenesslax                      0.000835 ***
## positionphrase-final:sylposN:tensenesslax 0.011938 *  
## positionphrase-final:sylposC:tensenesslax 0.108316    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- sylpsN sylpsC tnsnss pst-:N pst-:C pstn-: sylpN:
## pstnphrs-fn -0.018                                                        
## sylposN     -0.240  0.349                                                 
## sylposC     -0.240  0.349  0.500                                          
## tensenesslx -0.627 -0.122  0.192  0.192                                   
## pstnphrs-:N  0.171 -0.489 -0.713 -0.356 -0.137                            
## pstnphrs-:C  0.171 -0.489 -0.356 -0.713 -0.137  0.500                     
## pstnphrs-f: -0.136 -0.563 -0.311 -0.311  0.219  0.436  0.436              
## sylpsN:tnsn  0.171 -0.248 -0.710 -0.355 -0.270  0.506  0.253  0.438       
## sylpsC:tnsn  0.171 -0.248 -0.355 -0.710 -0.270  0.253  0.506  0.438  0.500
## pstnphr-:N: -0.121  0.347  0.505  0.253  0.192 -0.709 -0.354 -0.615 -0.712
## pstnphr-:C: -0.121  0.347  0.253  0.505  0.192 -0.354 -0.709 -0.615 -0.356
##             sylpC: ps-:N:
## pstnphrs-fn              
## sylposN                  
## sylposC                  
## tensenesslx              
## pstnphrs-:N              
## pstnphrs-:C              
## pstnphrs-f:              
## sylpsN:tnsn              
## sylpsC:tnsn              
## pstnphr-:N: -0.356       
## pstnphr-:C: -0.712  0.500
## convergence code: 0
## boundary (singular) fit: see ?isSingular
#models with additional random slopes do not converge

qqPlot(resid(m4))
## [1] 1589 1890
abline(h = cut<- 105)

mono2<-droplevels(subset(mono, resid(m4)<cut))
m4.1<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+position|subject) + (0+position|word), data = mono2)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 2 negative eigenvalues: -5.3e-03 -6.4e+01
# after cut model m4.1 does not converge any more

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


m5<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+position+tenseness|subject) + (0+position|word), data = mono2)
## boundary (singular) fit: see ?isSingular
m5.r2dummy<-lmer(sylposdur ~ position.num * sylpos.num * tenseness.num +
           (1+position.num+tenseness.num|subject) + (0+position.num|word), data = mono2)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0120507 (tol = 0.002, component 1)
r.squaredGLMM(m5.r2dummy)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.2914292 0.4106092
qqPlot(resid(m5))

## [1] 2085 1895
summary(m5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 |  
##     word) + (0 + position + tenseness | subject) + (0 + position |      word)
##    Data: mono2
## 
## REML criterion at convergence: 22891.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7908 -0.6844 -0.0448  0.5635  3.4316 
## 
## Random effects:
##  Groups    Name                  Variance  Std.Dev.  Corr       
##  subject   positionphrase-medial 0.000e+00  0.000000            
##            positionphrase-final  6.994e+01  8.363073   NaN      
##            tensenesslax          2.130e+00  1.459623   NaN -1.00
##  subject.1 (Intercept)           6.743e+01  8.211571            
##  word      positionphrase-medial 1.091e+02 10.442869            
##            positionphrase-final  2.084e+02 14.434746 1.00       
##  word.1    (Intercept)           7.259e-06  0.002694            
##  Residual                        1.025e+03 32.018253            
## Number of obs: 2342, groups:  subject, 10; word, 8
## 
## Fixed effects:
##                                           Estimate Std. Error       df t value
## (Intercept)                                 93.819      6.282   10.483  14.934
## positionphrase-final                        10.561      4.658   29.651   2.267
## sylposN                                      6.319      3.286 2304.606   1.923
## sylposC                                    -29.373      3.286 2304.606  -8.940
## tensenesslax                                -5.924      8.092    7.651  -0.732
## positionphrase-final:sylposN                34.143      4.608 2304.586   7.410
## positionphrase-final:sylposC                53.155      4.608 2304.587  11.536
## positionphrase-final:tensenesslax           -2.728      5.399   28.616  -0.505
## sylposN:tensenesslax                       -31.468      4.619 2304.477  -6.813
## sylposC:tensenesslax                        13.216      4.622 2304.498   2.859
## positionphrase-final:sylposN:tensenesslax  -16.057      6.488 2304.481  -2.475
## positionphrase-final:sylposC:tensenesslax   11.868      6.491 2304.491   1.829
##                                           Pr(>|t|)    
## (Intercept)                               2.11e-08 ***
## positionphrase-final                       0.03083 *  
## sylposN                                    0.05456 .  
## sylposC                                    < 2e-16 ***
## tensenesslax                               0.48593    
## positionphrase-final:sylposN              1.76e-13 ***
## positionphrase-final:sylposC               < 2e-16 ***
## positionphrase-final:tensenesslax          0.61727    
## sylposN:tensenesslax                      1.22e-11 ***
## sylposC:tensenesslax                       0.00428 ** 
## positionphrase-final:sylposN:tensenesslax  0.01340 *  
## positionphrase-final:sylposC:tensenesslax  0.06760 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- sylpsN sylpsC tnsnss pst-:N pst-:C pstn-: sylpN:
## pstnphrs-fn  0.170                                                        
## sylposN     -0.264  0.356                                                 
## sylposC     -0.264  0.356  0.505                                          
## tensenesslx -0.644 -0.164  0.205  0.205                                   
## pstnphrs-:N  0.188 -0.499 -0.713 -0.360 -0.146                            
## pstnphrs-:C  0.188 -0.499 -0.360 -0.713 -0.146  0.505                     
## pstnphrs-f: -0.146 -0.585 -0.307 -0.307  0.231  0.431  0.431              
## sylpsN:tnsn  0.188 -0.253 -0.711 -0.359 -0.287  0.507  0.256  0.431       
## sylpsC:tnsn  0.188 -0.253 -0.359 -0.711 -0.287  0.256  0.507  0.430  0.503
## pstnphr-:N: -0.134  0.354  0.506  0.256  0.204 -0.710 -0.358 -0.604 -0.712
## pstnphr-:C: -0.134  0.354  0.256  0.506  0.204 -0.358 -0.710 -0.604 -0.358
##             sylpC: ps-:N:
## pstnphrs-fn              
## sylposN                  
## sylposC                  
## tensenesslx              
## pstnphrs-:N              
## pstnphrs-:C              
## pstnphrs-f:              
## sylpsN:tnsn              
## sylpsC:tnsn              
## pstnphr-:N: -0.358       
## pstnphr-:C: -0.712  0.503
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# best model is m5


# model withouth stiel
no5<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+position+tenseness|subject) + (0+position|word), data = nostiel)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0103595 (tol = 0.002, component 1)
plot(allEffects(no5))

plot(allEffects(m5))

summary(no5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 |  
##     word) + (0 + position + tenseness | subject) + (0 + position |      word)
##    Data: nostiel
## 
## REML criterion at convergence: 16191.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7503 -0.6303 -0.0323  0.6253  4.1803 
## 
## Random effects:
##  Groups    Name                  Variance  Std.Dev. Corr       
##  subject   positionphrase-medial 1.167e+02 10.80416            
##            positionphrase-final  7.859e+01  8.86530  0.45      
##            tensenesslax          3.377e+00  1.83767 -1.00 -0.48
##  subject.1 (Intercept)           1.558e-04  0.01248            
##  word      positionphrase-medial 2.856e+01  5.34437            
##            positionphrase-final  4.081e+01  6.38845 1.00       
##  word.1    (Intercept)           3.834e-03  0.06192            
##  Residual                        5.557e+02 23.57381            
## Number of obs: 1767, groups:  subject, 10; word, 6
## 
## Fixed effects:
##                                            Estimate Std. Error        df
## (Intercept)                                 69.6699     5.0058   13.4724
## positionphrase-final                         8.2600     4.3434   20.9189
## sylposN                                     36.2687     2.7782 1732.9493
## sylposC                                     -2.6334     2.7782 1732.9493
## tensenesslax                                 0.1611     5.2041    6.3374
## positionphrase-final:sylposN                40.9826     3.8959 1732.9493
## positionphrase-final:sylposC                48.5343     3.8959 1732.9493
## positionphrase-final:tensenesslax           -4.4745     3.9804  116.3282
## sylposN:tensenesslax                       -41.2723     3.9222 1732.9493
## sylposC:tensenesslax                         8.4041     3.9222 1732.9493
## positionphrase-final:sylposN:tensenesslax  -18.8357     5.4958 1732.9493
## positionphrase-final:sylposC:tensenesslax   16.6428     5.4958 1732.9493
##                                           t value Pr(>|t|)    
## (Intercept)                                13.918 2.22e-09 ***
## positionphrase-final                        1.902 0.071065 .  
## sylposN                                    13.055  < 2e-16 ***
## sylposC                                    -0.948 0.343324    
## tensenesslax                                0.031 0.976262    
## positionphrase-final:sylposN               10.520  < 2e-16 ***
## positionphrase-final:sylposC               12.458  < 2e-16 ***
## positionphrase-final:tensenesslax          -1.124 0.263270    
## sylposN:tensenesslax                      -10.523  < 2e-16 ***
## sylposC:tensenesslax                        2.143 0.032276 *  
## positionphrase-final:sylposN:tensenesslax  -3.427 0.000624 ***
## positionphrase-final:sylposC:tensenesslax   3.028 0.002496 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- sylpsN sylpsC tnsnss pst-:N pst-:C pstn-: sylpN:
## pstnphrs-fn -0.431                                                        
## sylposN     -0.277  0.320                                                 
## sylposC     -0.277  0.320  0.500                                          
## tensenesslx -0.590  0.142  0.267  0.267                                   
## pstnphrs-:N  0.198 -0.448 -0.713 -0.357 -0.190                            
## pstnphrs-:C  0.198 -0.448 -0.357 -0.713 -0.190  0.500                     
## pstnphrs-f:  0.101 -0.460 -0.349 -0.349 -0.192  0.489  0.489              
## sylpsN:tnsn  0.197 -0.227 -0.708 -0.354 -0.377  0.505  0.253  0.493       
## sylpsC:tnsn  0.197 -0.227 -0.354 -0.708 -0.377  0.253  0.505  0.493  0.500
## pstnphr-:N: -0.140  0.318  0.506  0.253  0.269 -0.709 -0.354 -0.690 -0.714
## pstnphr-:C: -0.140  0.318  0.253  0.506  0.269 -0.354 -0.709 -0.690 -0.357
##             sylpC: ps-:N:
## pstnphrs-fn              
## sylposN                  
## sylposC                  
## tensenesslx              
## pstnphrs-:N              
## pstnphrs-:C              
## pstnphrs-f:              
## sylpsN:tnsn              
## sylpsC:tnsn              
## pstnphr-:N: -0.357       
## pstnphr-:C: -0.714  0.500
## convergence code: 0
## Model failed to converge with max|grad| = 0.0103595 (tol = 0.002, component 1)

1.3 Overall lengthening for each syllable constituent

# Without tenseness group
mono %>% group_by(sylpos, position) %>%
  summarise(mean.ms = mean(sylposdur)) %>%
  spread(position, mean.ms) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))
## `summarise()` regrouping output by 'sylpos' (override with `.groups` argument)
## # A tibble: 3 x 5
## # Groups:   sylpos [3]
##   sylpos medial final percent lengthening
##   <fct>   <dbl> <dbl>   <dbl> <chr>      
## 1 O        92.3  101.    1.10 10 %       
## 2 N        81.2  117.    1.44 44 %       
## 3 C        68.3  136.    2.00 100 %
# With tenseness group
mono %>% group_by(sylpos, position,tenseness) %>%
  summarise(mean.ms = mean(sylposdur)) %>%
  spread(position, mean.ms) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))
## `summarise()` regrouping output by 'sylpos', 'position' (override with `.groups` argument)
## # A tibble: 6 x 6
## # Groups:   sylpos [3]
##   sylpos tenseness medial final percent lengthening
##   <fct>  <fct>      <dbl> <dbl>   <dbl> <chr>      
## 1 O      tense       96.2 106.     1.11 11 %       
## 2 O      lax         88.4  96.2    1.09 9 %        
## 3 N      tense       99.9 145.     1.45 45 %       
## 4 N      lax         62.7  88.6    1.41 41 %       
## 5 C      tense       64.2 128.     2.00 100 %      
## 6 C      lax         72.3 144.     2.00 100 %
# normalized
mono %>% group_by(sylpos, position) %>%
  summarise(mean.norm = mean(sylposdurnorm)) %>%
  spread(position, mean.norm) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))
## `summarise()` regrouping output by 'sylpos' (override with `.groups` argument)
## # A tibble: 3 x 5
## # Groups:   sylpos [3]
##   sylpos medial final percent lengthening
##   <fct>   <dbl> <dbl>   <dbl> <chr>      
## 1 O       0.498 0.548    1.10 10 %       
## 2 N       0.435 0.627    1.44 44 %       
## 3 C       0.368 0.738    2.01 101 %
# Disyllabic without tenseness
di %>%  mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position, sylpos) %>%
  summarise(mean.ms = mean(sylposdur)) %>%
  spread(position, mean.ms) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))
## `summarise()` regrouping output by 'position' (override with `.groups` argument)
## # A tibble: 4 x 5
##   sylpos medial final percent lengthening
##   <fct>   <dbl> <dbl>   <dbl> <chr>      
## 1 O1       72.6  78.2    1.08 8 %        
## 2 N1       59.5  73.3    1.23 23 %       
## 3 O2       80.1  98.4    1.23 23 %       
## 4 N2       48.7 149.     3.06 206 %
# Disyllabic with tenseness
di %>%  mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position, sylpos,tenseness) %>%
  summarise(mean.ms = mean(sylposdur)) %>%
  spread(position, mean.ms) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))
## `summarise()` regrouping output by 'position', 'sylpos' (override with `.groups` argument)
## # A tibble: 8 x 6
## # Groups:   sylpos [4]
##   sylpos tenseness medial final percent lengthening
##   <fct>  <fct>      <dbl> <dbl>   <dbl> <chr>      
## 1 O1     tense       77.6  80.1    1.03 3 %        
## 2 O1     lax         67.6  76.3    1.13 13 %       
## 3 N1     tense       74.2  93.1    1.25 25 %       
## 4 N1     lax         44.9  53.3    1.19 19 %       
## 5 O2     tense       79.7  95.9    1.20 20 %       
## 6 O2     lax         80.4 101.     1.26 26 %       
## 7 N2     tense       44.0 147.     3.34 234 %      
## 8 N2     lax         53.2 151.     2.83 183 %

1.3.1 Excursus: Normalization

Normalization of dependent variable for speech rate does not change the effects considerably

m5.norm<-lmer(sylposdurnorm ~ position * sylpos * tenseness + (1|subject) + (1|word) +
                (0+position+sylpos|subject) + (0+position|word), data = mono)
## boundary (singular) fit: see ?isSingular
plot(allEffects(m5.norm))

plot(allEffects(m5))

1.3.2 Figure 2 (normalized stacked mean duration per tenseness and condition)

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

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

Figure 3: Acoustic mean duration of syllable constituents per phrasal position

# Duration in milliseconds
(plot.ac.mono<- mono %>% group_by(position,word,tenseness, sylpos) %>%
  dplyr::summarise(mean.dur = mean(sylposdur), sd.dur = sd(sylposdur)) %>%
  filter(position=="phrase-final") %>%
  ungroup() %>%
  mutate(lengthening = leng.perc$lengthening) %>%
  group_by(word) %>%
  mutate(ymin=cumsum(mean.dur) -sd.dur/2, ymax=cumsum(mean.dur) + sd.dur/2) %>%
  ggplot(aes(x=word, y=mean.dur, fill=sylpos))+
  geom_bar(
    stat="identity",
    position=position_stack(reverse=TRUE),alpha=0.9
  )+
  geom_errorbar(aes(ymin=ymin, ymax=ymax),
                stat="identity",
                width=.2                    # Width of the error bars
  )+
  coord_flip()+
  ylab("")+
  geom_text(aes(x=word, y=mean.dur, label = lengthening),
            color="black", hjust= 2,size = 3, position=position_stack(reverse=TRUE))+
  annotate("text", x=1:8,y= 500,label=rep(c("tense","lax"),4),size=3)+
  # ylim(0,2.7)+
  # scale_fill_grey(start=0,end=0.6, name="")+
  xlab("Monosyllabic targets")+
  scale_fill_brewer(palette="Set3",name="",labels=c("Onset","Nucleus","Coda"))+
  theme_bw(base_size = 10) +
  theme( legend.position="top"))
## `summarise()` regrouping output by 'position', 'word', 'tenseness' (override with `.groups` argument)

(plot.ac.di<-
  di %>% mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position,word,tenseness, sylpos) %>%
  summarise(mean.dur = mean(sylposdur), sd.dur = sd(sylposdur)) %>%
  filter(position=="phrase-final") %>%
  ungroup() %>%
  mutate(lengthening = leng.perc.di$lengthening) %>%
  group_by(word) %>%
  mutate(ymin=cumsum(mean.dur) -sd.dur/2, ymax=cumsum(mean.dur) + sd.dur/2) %>%
  ggplot(aes(x=word, y=mean.dur, fill=sylpos))+
  geom_bar(
    stat="identity",
    position=position_stack(reverse=T),alpha=0.9
  )+
  geom_errorbar(aes(ymin=ymin, ymax=ymax),
                stat="identity",
                width=.2                    # Width of the error bars
  )+
  coord_flip()+
  ylab("Acoustic duration (ms)")+
  geom_text(aes(x=word, y=mean.dur, label = lengthening),
            color="black", vjust= 0,hjust=2,size = 3, position=position_stack(reverse=T))+
  annotate("text", x=1:4,y= 500,label=rep(c("tense","lax"),2),size=3)+
  # scale_fill_grey(start=0,end=0.6, name="")+
  xlab("Disyllabic targets")+
  scale_fill_brewer(palette="Set2",name="",labels=c("Onset 1","Nucleus 1","Onset 2","Nucleus 2"))+
  theme_bw(base_size = 10) +
  theme( legend.position="top"))
## `summarise()` regrouping output by 'position', 'word', 'tenseness' (override with `.groups` argument)

ggpubr::ggarrange(plot.ac.mono, plot.ac.di, common.legend = F,labels = c("a)","b)"),
          font.label = list(size=11,face="plain"),align = "h", widths = c(1.1,1), ncol = 1)

ggsave(file="v2.stackedacousticfinal.pdf", width = 200, height = 160, units = "mm")


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

(plot.ac.mono.data<-
  ggplot(mono %>% filter(sylpos=="N"), aes(x=tenseness, y=sylposdur,fill=position))+
  geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),alpha=.3,size=.7,col="black")+
  geom_boxplot(alpha=.3)+facet_wrap(~wordpair,nrow=1)+
  scale_fill_brewer(palette = "Set1")+
  xlab("Tenseness")+
  ylab("Acoustic duration (ms)")+
  theme_bw(base_size=10)+
  theme(legend.position = "top",legend.title = element_blank())
)

(plot.ac.di.data<-ggplot(di %>% filter(sylpos=="N1"), aes(x=tenseness, y=sylposdur, fill=position))+
  geom_boxplot(alpha=.3)+
  geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),alpha=.3,size=.7,col="black")+
  facet_wrap(~wordpair,nrow=1)+
  scale_fill_brewer(palette = "Set1")+
  xlab("Tenseness")+
  ylab("")+
  theme_bw(base_size=10)+
  theme(legend.position = "top",legend.title = element_blank())
)

cairo_pdf(filename = "v2.ac-vowel-tenseness-position.pdf",width=8,height=5)
ggpubr::ggarrange(plot.ac.mono.data, plot.ac.di.data, common.legend = T,labels = c("a)","b)"),
                  font.label = list(size=11,face="plain"),align = "h", widths = c(1.1,1))
dev.off()
## png 
##   2

1.3.3 Subset model acoustic monosyllabic for lax only

mono.lax<-droplevels(subset(mono, tenseness=="lax"))
mono.lax.nucl<-droplevels(subset(mono, tenseness=="lax" & sylpos=="N"))
mono.tense<-droplevels(subset(mono, tenseness=="tense"))
mono.tense.nucl<-droplevels(subset(mono, tenseness=="tense" & sylpos=="N"))

l1<-lmer(sylposdur ~ position + sylpos  + (1|subject) + (1|word), data = mono.lax)
l2<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word), data = mono.lax)
anova(l1,l2)
## refitting model(s) with ML (instead of REML)
## Data: mono.lax
## Models:
## l1: sylposdur ~ position + sylpos + (1 | subject) + (1 | word)
## l2: sylposdur ~ position * sylpos + (1 | subject) + (1 | word)
##    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## l1  7 11590 11626 -5788.1    11576                             
## l2  9 11359 11405 -5670.4    11341 235.22      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
l3<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax)
## boundary (singular) fit: see ?isSingular
anova(l2,l3,refit=F) #better
## Data: mono.lax
## Models:
## l2: sylposdur ~ position * sylpos + (1 | subject) + (1 | word)
## l3: sylposdur ~ position * sylpos + (1 | subject) + (1 | word) + 
## l3:     (0 + position + sylpos | word)
##    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## l2  9 11335 11381 -5658.7    11317                             
## l3 19 10448 10544 -5204.8    10410 907.61     10  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(resid(l3))
## [1] 844 530
abline(h = cut <- 50)

mono.lax.sub<-droplevels(subset(mono.lax, resid(l3) < cut))
l4<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax.sub)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(l4))
## [1] 1155  692
abline(h = cut <- -50)

mono.lax.sub<-droplevels(subset(mono.lax.sub, resid(l4) > cut))
l5<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax.sub)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(l5))

## [1] 560 797
plot(allEffects(l5))

summary(l5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position * sylpos + (1 | subject) + (1 | word) +  
##     (0 + position + sylpos | word)
##    Data: mono.lax.sub
## 
## REML criterion at convergence: 9859.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.89003 -0.65202 -0.01035  0.59712  3.05073 
## 
## Random effects:
##  Groups   Name                  Variance  Std.Dev.  Corr             
##  subject  (Intercept)           5.174e+01 7.193e+00                  
##  word     (Intercept)           3.886e-07 6.234e-04                  
##  word.1   positionphrase-medial 1.741e+03 4.172e+01                  
##           positionphrase-final  2.228e+03 4.720e+01  1.00            
##           sylposN               2.795e+03 5.287e+01 -0.98 -0.98      
##           sylposC               2.400e+03 4.899e+01 -1.00 -1.00  0.98
##  Residual                       2.872e+02 1.695e+01                  
## Number of obs: 1155, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                    86.896     21.020    3.062   4.134   0.0247 *  
## positionphrase-final            9.554      3.343    4.429   2.858   0.0408 *  
## sylposN                       -23.989     26.492    2.997  -0.906   0.4320    
## sylposC                       -14.912     24.555    2.999  -0.607   0.5866    
## positionphrase-final:sylposN   15.762      2.425 1131.024   6.500  1.2e-10 ***
## positionphrase-final:sylposC   60.525      2.459 1131.316  24.612  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- sylpsN sylpsC pst-:N
## pstnphrs-fn  0.787                            
## sylposN     -0.972 -0.790                     
## sylposC     -0.991 -0.777  0.974              
## pstnphrs-:N  0.029 -0.365 -0.046 -0.025       
## pstnphrs-:C  0.029 -0.360 -0.023 -0.049  0.496
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lattice::dotplot(ranef(l5))
## $subject

## 
## $word

1.3.4 Proportional lengthening of vowels in nuclei, monosyllabic, with and without Stiel/still

tn1<-lmer(sylposdur ~ position  + (1|subject) + (0+position|word), data = mono.tense.nucl)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0297967 (tol = 0.002, component 1)
qqPlot(resid(tn1))

## [1] 254 278
summary(tn1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position + (1 | subject) + (0 + position | word)
##    Data: mono.tense.nucl
## 
## REML criterion at convergence: 3464.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5604 -0.6538 -0.1107  0.6794  3.1496 
## 
## Random effects:
##  Groups   Name                  Variance Std.Dev. Corr
##  subject  (Intercept)            315.1   17.75        
##  word     positionphrase-medial  664.3   25.77        
##           positionphrase-final  1190.7   34.51    0.99
##  Residual                        379.0   19.47        
## Number of obs: 390, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                      Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)           100.508     14.127   4.191   7.115  0.00173 **
## positionphrase-final   43.895      5.043   3.011   8.704  0.00315 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pstnphrs-fn 0.745 
## convergence code: 0
## Model failed to converge with max|grad| = 0.0297967 (tol = 0.002, component 1)
tn1.coef<-summary(tn1)$coefficients
# Proportion lengthening of tense vowel with Stiel
(tn1.coef[1,1]+tn1.coef[2,1])/tn1.coef[1,1] # tense with Stiel
## [1] 1.436735
# Model of proportional lengthening in monosyllabic nuclei
mono.proplen <- mono %>%
  filter(sylpos=="N") %>%
  group_by(subject,tenseness,wordpair,position) %>%
  summarise(mean = mean(sylposdur)) %>%
  tidyr::spread(position, mean) %>%
  mutate(proplen = ((`phrase-final`/`phrase-medial`)-1)*100)
## `summarise()` regrouping output by 'subject', 'tenseness', 'wordpair' (override with `.groups` argument)
mp1 <- lmer(proplen ~ tenseness + (1|subject) + (1|wordpair), data=mono.proplen)
qqPlot(resid(mp1))

## [1] 80 52
summary(mp1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: proplen ~ tenseness + (1 | subject) + (1 | wordpair)
##    Data: mono.proplen
## 
## REML criterion at convergence: 722.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.85764 -0.55990 -0.06001  0.50219  2.36290 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 151.05   12.290  
##  wordpair (Intercept)  37.99    6.164  
##  Residual             466.85   21.607  
## Number of obs: 80, groups:  subject, 10; wordpair, 4
## 
## Fixed effects:
##              Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)    48.268      6.023 10.353   8.014 9.37e-06 ***
## tensenesslax   -5.756      4.831 66.000  -1.191    0.238    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tensenesslx -0.401
plot(allEffects(mp1))

# lax only
ln1<-lmer(sylposdur ~ position  + (1|subject) + (1|word), data = mono.lax.nucl)
qqPlot(resid(ln1))
## [1] 328  33
abline(h = cut <- c(-40,40))

mono.lax.nucl.sub<-droplevels(subset(mono.lax.nucl, resid(ln1) > cut[1] & resid(ln1) < cut[2]))
ln2<-lmer(sylposdur ~ position  + (1|subject) + (1|word), data = mono.lax.nucl.sub)
qqPlot(resid(ln2))

## [1] 387  36
summary(ln2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position + (1 | subject) + (1 | word)
##    Data: mono.lax.nucl.sub
## 
## REML criterion at convergence: 3063.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1331 -0.5632 -0.0549  0.5202  3.3720 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 127.1    11.27   
##  word     (Intercept) 160.4    12.66   
##  Residual             138.4    11.76   
## Number of obs: 389, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                      Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)            62.931      7.315   5.075   8.603 0.000325 ***
## positionphrase-final   25.759      1.195 375.053  21.559  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pstnphrs-fn -0.082
ln2.coef<-summary(ln2)$coefficients
# Proportion lengthening of lax vowels with Stiel
(ln2.coef[1,1]+ln2.coef[2,1])/ln2.coef[1,1]
## [1] 1.409323
#without Stiel/stiel tense
mono.tense.nucl.nostiel<-droplevels(subset(mono, tenseness=="tense" & sylpos=="N" & !(word%in%c("Stiel","still")) ))
mono.lax.nucl.nostiel<-droplevels(subset(mono, tenseness=="lax" & sylpos=="N" & !(word%in%c("Stiel","still")) ))

#tense
tn.nostiel<-lmer(sylposdur ~ position  + (1|subject) + (0+position|word), data = mono.tense.nucl.nostiel)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(tn.nostiel))

## [1] 196 208
summary(tn.nostiel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position + (1 | subject) + (0 + position | word)
##    Data: mono.tense.nucl.nostiel
## 
## REML criterion at convergence: 2592.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6745 -0.6306 -0.1063  0.6545  3.1057 
## 
## Random effects:
##  Groups   Name                  Variance Std.Dev. Corr
##  subject  (Intercept)            367.6   19.17        
##  word     positionphrase-medial  798.6   28.26        
##           positionphrase-final  1154.6   33.98    1.00
##  Residual                        362.8   19.05        
## Number of obs: 293, groups:  subject, 10; word, 3
## 
## Fixed effects:
##                      Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)           106.437     17.478   2.573    6.09  0.01366 * 
## positionphrase-final   48.310      3.984   2.211   12.13  0.00459 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pstnphrs-fn 0.738 
## convergence code: 0
## boundary (singular) fit: see ?isSingular
tn1.coef.nostiel<-summary(tn.nostiel)$coefficients
# Proportion lengthening of tense vowel without Stiel
(tn1.coef.nostiel[1,1]+tn1.coef.nostiel[2,1])/tn1.coef.nostiel[1,1]
## [1] 1.45389
#lax
ln.nostiel<-lmer(sylposdur ~ position  + (1|subject) + (0+position|word), data = mono.lax.nucl.nostiel)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(ln.nostiel))

## [1] 254 220
summary(ln.nostiel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position + (1 | subject) + (0 + position | word)
##    Data: mono.lax.nucl.nostiel
## 
## REML criterion at convergence: 2276.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0992 -0.5366 -0.0712  0.5662  2.9576 
## 
## Random effects:
##  Groups   Name                  Variance Std.Dev. Corr
##  subject  (Intercept)           118.0    10.86        
##  word     positionphrase-medial 229.1    15.14        
##           positionphrase-final  202.7    14.24    1.00
##  Residual                       113.2    10.64        
## Number of obs: 296, groups:  subject, 10; word, 3
## 
## Fixed effects:
##                      Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)            65.258      9.432  2.644   6.919   0.0092 ** 
## positionphrase-final   25.298      1.344  7.053  18.816 2.74e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pstnphrs-fn -0.420
## convergence code: 0
## boundary (singular) fit: see ?isSingular
ln1.coef.nostiel<-summary(ln.nostiel)$coefficients
# Proportion lengthening of lax vowel without Stiel
(ln1.coef.nostiel[1,1]+ln1.coef.nostiel[2,1])/ln1.coef.nostiel[1,1]
## [1] 1.387658
l2<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word), data = mono.lax)
anova(l1,l2)
## refitting model(s) with ML (instead of REML)
## Data: mono.lax
## Models:
## l1: sylposdur ~ position + sylpos + (1 | subject) + (1 | word)
## l2: sylposdur ~ position * sylpos + (1 | subject) + (1 | word)
##    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## l1  7 11590 11626 -5788.1    11576                             
## l2  9 11359 11405 -5670.4    11341 235.22      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
l3<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax)
## boundary (singular) fit: see ?isSingular
anova(l2,l3,refit=F) #better
## Data: mono.lax
## Models:
## l2: sylposdur ~ position * sylpos + (1 | subject) + (1 | word)
## l3: sylposdur ~ position * sylpos + (1 | subject) + (1 | word) + 
## l3:     (0 + position + sylpos | word)
##    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## l2  9 11335 11381 -5658.7    11317                             
## l3 19 10448 10544 -5204.8    10410 907.61     10  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(resid(l3))
## [1] 844 530
abline(h = cut <- 50)

mono.lax.sub<-droplevels(subset(mono.lax, resid(l3) < cut))
l4<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax.sub)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(l4))
## [1] 1155  692
abline(h = cut <- -50)

mono.lax.sub<-droplevels(subset(mono.lax.sub, resid(l4) > cut))
l5<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax.sub)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(l5))

## [1] 560 797
plot(allEffects(l5))

summary(l5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position * sylpos + (1 | subject) + (1 | word) +  
##     (0 + position + sylpos | word)
##    Data: mono.lax.sub
## 
## REML criterion at convergence: 9859.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.89003 -0.65202 -0.01035  0.59712  3.05073 
## 
## Random effects:
##  Groups   Name                  Variance  Std.Dev.  Corr             
##  subject  (Intercept)           5.174e+01 7.193e+00                  
##  word     (Intercept)           3.886e-07 6.234e-04                  
##  word.1   positionphrase-medial 1.741e+03 4.172e+01                  
##           positionphrase-final  2.228e+03 4.720e+01  1.00            
##           sylposN               2.795e+03 5.287e+01 -0.98 -0.98      
##           sylposC               2.400e+03 4.899e+01 -1.00 -1.00  0.98
##  Residual                       2.872e+02 1.695e+01                  
## Number of obs: 1155, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                    86.896     21.020    3.062   4.134   0.0247 *  
## positionphrase-final            9.554      3.343    4.429   2.858   0.0408 *  
## sylposN                       -23.989     26.492    2.997  -0.906   0.4320    
## sylposC                       -14.912     24.555    2.999  -0.607   0.5866    
## positionphrase-final:sylposN   15.762      2.425 1131.024   6.500  1.2e-10 ***
## positionphrase-final:sylposC   60.525      2.459 1131.316  24.612  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- sylpsN sylpsC pst-:N
## pstnphrs-fn  0.787                            
## sylposN     -0.972 -0.790                     
## sylposC     -0.991 -0.777  0.974              
## pstnphrs-:N  0.029 -0.365 -0.046 -0.025       
## pstnphrs-:C  0.029 -0.360 -0.023 -0.049  0.496
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lattice::dotplot(ranef(l5))
## $subject

## 
## $word

# Model of proportional lengthening in disyllabic nuclei
di.proplen <- di %>%
  filter(sylpos=="N1") %>%
  group_by(subject,tenseness,wordpair,position) %>%
  summarise(mean = mean(sylposdur)) %>%
  tidyr::spread(position, mean) %>%
  mutate(proplen = ((`phrase-final`/`phrase-medial`)-1)*100)
## `summarise()` regrouping output by 'subject', 'tenseness', 'wordpair' (override with `.groups` argument)
dp1 <- lmer(proplen ~ tenseness + (1|subject) + (1|wordpair), data=di.proplen)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(dp1))

## [1] 21 32
summary(dp1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: proplen ~ tenseness + (1 | subject) + (1 | wordpair)
##    Data: di.proplen
## 
## REML criterion at convergence: 351.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3945 -0.7795 -0.1186  0.6708  3.0899 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept)   0.0     0.00   
##  wordpair (Intercept)   0.0     0.00   
##  Residual             521.8    22.84   
## Number of obs: 40, groups:  subject, 10; wordpair, 2
## 
## Fixed effects:
##              Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)    29.847      5.108 38.000   5.843 9.34e-07 ***
## tensenesslax   -6.635      7.224 38.000  -0.918    0.364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tensenesslx -0.707
## convergence code: 0
## boundary (singular) fit: see ?isSingular
plot(allEffects(dp1))

# Table with proplen
x<-capture.output(
  texreg::texreg(list(texreg::extract(mp1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                      texreg::extract(dp1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                      # texreg::extract(mf3.1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
  ),
  digits = 1,booktabs = T,custom.model.names = c("Monosyllabic \\%","Disyllabic \\%"),
  float.pos = "!htbp", label = "tab:proplen", use.packages = F,single.row = T,caption.above = T,
  caption="Model predictions for the proportional lengthening of vowels  in  mono- and disyllabic words
   (standards error in parentheses) for the fixed effects tenseness (tense vs. lax) and random intercepts for subjects and word pairs.",
  dcolumn = T, fontsize = "scriptsize",center=T)
)
x[5]<-paste0(x[5], "\\centering")
x[16]<-paste0(x[16],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(mp1)[1],2),"/",round(r.squaredGLMM(mp1)[2],2),
              "&", round(r.squaredGLMM(dp1)[1],2), "/", round(r.squaredGLMM(dp1)[2],2),
              "\\\\ \n")
cat(x, sep = "\n")
## 
## \begin{table}[!htbp]
## \caption{Model predictions for the proportional lengthening of vowels  in  mono- and disyllabic words
##    (standards error in parentheses) for the fixed effects tenseness (tense vs. lax) and random intercepts for subjects and word pairs.}
## \begin{center}\centering
## \begin{scriptsize}
## \begin{tabular}{l D{)}{)}{9)3} D{)}{)}{9)3} }
## \toprule
##  & \multicolumn{1}{c}{Monosyllabic \%} & \multicolumn{1}{c}{Disyllabic \%} \\
## \midrule
## (Intercept)           & 48.3 \; (6.0)^{***} & 29.8 \; (5.1)^{***} \\
## tensenesslax          & -5.8 \; (4.8)       & -6.6 \; (7.2)       \\
## \midrule
## AIC                   & 732.5               & 361.6               \\
## Num. obs.             & 80                  & 40                  \\
## Num. groups: subject  & 10                  & 10                  \\\midrule R$_m^2$/R$_c^2$ & 0.01/0.3&0.02/0.02\\ 
## 
## Num. groups: wordpair & 4                   & 2                   \\
## \bottomrule
## \multicolumn{3}{l}{\tiny{$^{***}p<0.001$, $^{**}p<0.01$, $^*p<0.05$}}
## \end{tabular}
## \end{scriptsize}
## \label{tab:proplen}
## \end{center}
## \end{table}

1.3.5 Post-hoc tests (multiple comparisons, Tukey-corrected) for mono, lax

l5em <- emmeans(l5, pairwise~position*sylpos)
l5em.tab<-l5em$contrasts %>%
  summary(infer = T)
l5em.tab$p.value<-lapply(l5em.tab$p.value, FUN = function (x) {
  # levels of p-values
  if (x>=0.05)  p <- "n.\\,s." # not significant
  if (x<0.05)  p <- "<0,05"  # significant
  if (x<=0.01)  p <- "<0,01"  # significant
  if (x<=0.001) p <- "<0,001" # significant
  return(p)
})
l5em.tab<-as.data.frame(l5em.tab)
colnames(l5em.tab)<-c("Comparison","$\\beta$","s","df","Lower CI","Upper CI","t","p")
l5em.tab.sig<-subset(l5em.tab, p!="n.\\,s.")
l5em.tab.sig
##                          Comparison  $\\beta$        s       df   Lower CI
## 10 phrase-medial,N - phrase-final,N -25.31610 3.338346 4.386355  -40.42224
## 12 phrase-medial,N - phrase-final,C -79.15585 7.598238 3.056582 -121.75405
## 14  phrase-final,N - phrase-final,C -53.83975 6.087517 3.267556  -86.42379
## 15 phrase-medial,C - phrase-final,C -70.07870 3.363707 4.517218  -85.08832
##     Upper CI          t      p
## 10 -10.20995  -7.583425  <0,01
## 12 -36.55764 -10.417659  <0,01
## 14 -21.25571  -8.844288  <0,05
## 15 -55.06908 -20.833769 <0,001

1.3.6 Monosyllabic Effect plots

Preparation

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

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

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

1.4 Disyllabic target words

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

d1<-lmer(sylposdur ~ position + sylpos + tenseness + (1|subject) + (1|word), data = di)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(d1))

## [1] 1402 1097
d2<-lmer(sylposdur ~ position * sylpos + tenseness + (1|subject) + (1|word), data = di)
## boundary (singular) fit: see ?isSingular
anova(d1,d2)
## refitting model(s) with ML (instead of REML)
## Data: di
## Models:
## d1: sylposdur ~ position + sylpos + tenseness + (1 | subject) + (1 | 
## d1:     word)
## d2: sylposdur ~ position * sylpos + tenseness + (1 | subject) + (1 | 
## d2:     word)
##    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## d1  9 15024 15072 -7502.9    15006                             
## d2 12 14123 14188 -7049.6    14099 906.56      3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d3<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word), data = di)
## boundary (singular) fit: see ?isSingular
anova(d2,d3)
## refitting model(s) with ML (instead of REML)
## Data: di
## Models:
## d2: sylposdur ~ position * sylpos + tenseness + (1 | subject) + (1 | 
## d2:     word)
## d3: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 | 
## d3:     word)
##    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## d2 12 14123 14188 -7049.6    14099                             
## d3 19 13889 13991 -6925.7    13851 247.84      7  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d4<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+tenseness|subject) + (0+position+tenseness+sylpos|word), data = di)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -2.0e+00
anova(d3,d4,refit=F)
## Data: di
## Models:
## d3: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 | 
## d3:     word)
## d4: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 | 
## d4:     word) + (0 + tenseness | subject) + (0 + position + tenseness + 
## d4:     sylpos | word)
##    Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## d3 19 13834 13936 -6898.1    13796                             
## d4 43 13681 13911 -6797.4    13595 201.36     24  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#best model is d4, models with additional random slopes do not converge
qqPlot(resid(d4))
## [1] 1402  731
abline(h = cut<- c(-50,60))

di2<-droplevels(subset(di, resid(d4)<cut[2] & resid(d4)> cut[1]))

di2$sylpos <- factor(di2$sylpos, levels = c("O1","N1","O2","N2"))
d5<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
             (0+tenseness|subject) + (0+position+tenseness+sylpos|word), data = di2)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.9e-01
d5.r2dummy<-lmer(sylposdur ~ position * sylpos * tenseness +
               (1+tenseness|subject) + (1+position+tenseness+sylpos|word), data = di2)
## boundary (singular) fit: see ?isSingular
anova(d5, d5.r2dummy, refit=F)
## Data: di2
## Models:
## d5.r2dummy: sylposdur ~ position * sylpos * tenseness + (1 + tenseness | 
## d5.r2dummy:     subject) + (1 + position + tenseness + sylpos | word)
## d5: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 | 
## d5:     word) + (0 + tenseness | subject) + (0 + position + tenseness + 
## d5:     sylpos | word)
##            Df   AIC   BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## d5.r2dummy 41 12854 13073 -6385.9    12772                        
## d5         43 12858 13088 -6385.9    12772     0      2          1
# after cut model does not converge any more
qqPlot(resid(d5))

## [1] 773 335
summary(d5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 |  
##     word) + (0 + tenseness | subject) + (0 + position + tenseness +  
##     sylpos | word)
##    Data: di2
## 
## REML criterion at convergence: 12771.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5127 -0.6642 -0.0445  0.5997  3.9051 
## 
## Random effects:
##  Groups    Name                  Variance  Std.Dev.  Corr                   
##  subject   tensenesstense        6.551e+01 8.094e+00                        
##            tensenesslax          3.947e+01 6.283e+00 0.96                   
##  subject.1 (Intercept)           3.400e+00 1.844e+00                        
##  word      positionphrase-medial 3.967e+02 1.992e+01                        
##            positionphrase-final  2.968e+02 1.723e+01  1.00                  
##            tensenesslax          1.775e-01 4.213e-01 -1.00 -1.00            
##            sylposN1              7.060e+02 2.657e+01 -1.00 -1.00  1.00      
##            sylposO2              3.043e+02 1.745e+01 -1.00 -1.00  1.00  1.00
##            sylposN2              6.696e+02 2.588e+01 -1.00 -1.00  1.00  1.00
##  word.1    (Intercept)           5.103e-07 7.144e-04                        
##  Residual                        2.346e+02 1.532e+01                        
##       
##       
##       
##       
##       
##       
##       
##       
##       
##   1.00
##       
##       
## Number of obs: 1541, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                                             Estimate Std. Error        df
## (Intercept)                                  75.7441    14.4127    1.8835
## positionphrase-final                          4.6753     2.9112    5.8053
## sylposN1                                     -1.6831    18.9173    1.7739
## sylposO2                                      3.7672    12.5325    1.8074
## sylposN2                                    -31.9499    18.4312    1.7758
## tensenesslax                                 -7.7560    19.8502    1.7714
## positionphrase-final:sylposN1                14.4965     3.1029 1502.9286
## positionphrase-final:sylposO2                11.8324     3.1028 1502.8700
## positionphrase-final:sylposN2                99.9974     3.1414 1503.6649
## positionphrase-final:tensenesslax             3.4240     4.1106    5.7696
## sylposN1:tensenesslax                       -21.3373    26.7512    1.7734
## sylposO2:tensenesslax                         8.7856    17.7208    1.8062
## sylposN2:tensenesslax                        15.2809    26.0637    1.7752
## positionphrase-final:sylposN1:tensenesslax  -14.1171     4.3823 1502.9774
## positionphrase-final:sylposO2:tensenesslax    0.4459     4.3820 1502.8573
## positionphrase-final:sylposN2:tensenesslax  -11.1218     4.4619 1505.2441
##                                            t value Pr(>|t|)    
## (Intercept)                                  5.255 0.038955 *  
## positionphrase-final                         1.606 0.161052    
## sylposN1                                    -0.089 0.938113    
## sylposO2                                     0.301 0.794730    
## sylposN2                                    -1.733 0.240638    
## tensenesslax                                -0.391 0.737913    
## positionphrase-final:sylposN1                4.672 3.25e-06 ***
## positionphrase-final:sylposO2                3.813 0.000143 ***
## positionphrase-final:sylposN2               31.832  < 2e-16 ***
## positionphrase-final:tensenesslax            0.833 0.437985    
## sylposN1:tensenesslax                       -0.798 0.517820    
## sylposO2:tensenesslax                        0.496 0.673740    
## sylposN2:tensenesslax                        0.586 0.623506    
## positionphrase-final:sylposN1:tensenesslax  -3.221 0.001303 ** 
## positionphrase-final:sylposO2:tensenesslax   0.102 0.918966    
## positionphrase-final:sylposN2:tensenesslax  -2.493 0.012788 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## convergence code: 0
## boundary (singular) fit: see ?isSingular

1.4.1 Excursus: Normalization

Normalization of dependent variable for speech rate does not change the effects considerably

d5.norm<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
                (0+tenseness|subject) + (0+position+tenseness+sylpos|word), data = di2)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.9e-01
plot(allEffects(d5.norm))

summary(d5.norm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position * sylpos * tenseness + (1 | subject) + (1 |  
##     word) + (0 + tenseness | subject) + (0 + position + tenseness +  
##     sylpos | word)
##    Data: di2
## 
## REML criterion at convergence: 12771.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5127 -0.6642 -0.0445  0.5997  3.9051 
## 
## Random effects:
##  Groups    Name                  Variance  Std.Dev.  Corr                   
##  subject   tensenesstense        6.551e+01 8.094e+00                        
##            tensenesslax          3.947e+01 6.283e+00 0.96                   
##  subject.1 (Intercept)           3.400e+00 1.844e+00                        
##  word      positionphrase-medial 3.967e+02 1.992e+01                        
##            positionphrase-final  2.968e+02 1.723e+01  1.00                  
##            tensenesslax          1.775e-01 4.213e-01 -1.00 -1.00            
##            sylposN1              7.060e+02 2.657e+01 -1.00 -1.00  1.00      
##            sylposO2              3.043e+02 1.745e+01 -1.00 -1.00  1.00  1.00
##            sylposN2              6.696e+02 2.588e+01 -1.00 -1.00  1.00  1.00
##  word.1    (Intercept)           5.103e-07 7.144e-04                        
##  Residual                        2.346e+02 1.532e+01                        
##       
##       
##       
##       
##       
##       
##       
##       
##       
##   1.00
##       
##       
## Number of obs: 1541, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                                             Estimate Std. Error        df
## (Intercept)                                  75.7441    14.4127    1.8835
## positionphrase-final                          4.6753     2.9112    5.8053
## sylposN1                                     -1.6831    18.9173    1.7739
## sylposO2                                      3.7672    12.5325    1.8074
## sylposN2                                    -31.9499    18.4312    1.7758
## tensenesslax                                 -7.7560    19.8502    1.7714
## positionphrase-final:sylposN1                14.4965     3.1029 1502.9286
## positionphrase-final:sylposO2                11.8324     3.1028 1502.8700
## positionphrase-final:sylposN2                99.9974     3.1414 1503.6649
## positionphrase-final:tensenesslax             3.4240     4.1106    5.7696
## sylposN1:tensenesslax                       -21.3373    26.7512    1.7734
## sylposO2:tensenesslax                         8.7856    17.7208    1.8062
## sylposN2:tensenesslax                        15.2809    26.0637    1.7752
## positionphrase-final:sylposN1:tensenesslax  -14.1171     4.3823 1502.9774
## positionphrase-final:sylposO2:tensenesslax    0.4459     4.3820 1502.8573
## positionphrase-final:sylposN2:tensenesslax  -11.1218     4.4619 1505.2441
##                                            t value Pr(>|t|)    
## (Intercept)                                  5.255 0.038955 *  
## positionphrase-final                         1.606 0.161052    
## sylposN1                                    -0.089 0.938113    
## sylposO2                                     0.301 0.794730    
## sylposN2                                    -1.733 0.240638    
## tensenesslax                                -0.391 0.737913    
## positionphrase-final:sylposN1                4.672 3.25e-06 ***
## positionphrase-final:sylposO2                3.813 0.000143 ***
## positionphrase-final:sylposN2               31.832  < 2e-16 ***
## positionphrase-final:tensenesslax            0.833 0.437985    
## sylposN1:tensenesslax                       -0.798 0.517820    
## sylposO2:tensenesslax                        0.496 0.673740    
## sylposN2:tensenesslax                        0.586 0.623506    
## positionphrase-final:sylposN1:tensenesslax  -3.221 0.001303 ** 
## positionphrase-final:sylposO2:tensenesslax   0.102 0.918966    
## positionphrase-final:sylposN2:tensenesslax  -2.493 0.012788 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## convergence code: 0
## boundary (singular) fit: see ?isSingular
###

1.4.2 Post-hoc tests (multiple comparisons, Tukey-corrected):

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

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

1.4.3 Subset model acoustic disyllabic for lax only

di.lax<-droplevels(subset(di, tenseness=="lax"))
dl1<-lmer(sylposdur ~ position + sylpos  + (1|subject) + (1|word), data = di.lax)
## boundary (singular) fit: see ?isSingular
dl2<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word), data = di.lax)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00259991 (tol = 0.002, component 1)
anova(dl1,dl2)
## refitting model(s) with ML (instead of REML)
## Data: di.lax
## Models:
## dl1: sylposdur ~ position + sylpos + (1 | subject) + (1 | word)
## dl2: sylposdur ~ position * sylpos + (1 | subject) + (1 | word)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## dl1  8 7448.1 7485.4 -3716.1   7432.1                             
## dl2 11 6968.6 7019.9 -3473.3   6946.6 485.49      3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dl3<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = di.lax)
## boundary (singular) fit: see ?isSingular
anova(dl2,dl3,refit=F) #better
## Data: di.lax
## Models:
## dl2: sylposdur ~ position * sylpos + (1 | subject) + (1 | word)
## dl3: sylposdur ~ position * sylpos + (1 | subject) + (1 | word) + 
## dl3:     (0 + position + sylpos | word)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## dl2 11 6940.1 6991.4 -3459.1   6918.1                             
## dl3 26 6895.5 7016.8 -3421.7   6843.5 74.614     15  6.649e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(resid(l3))
## [1] 844 530
abline(h = cut <- 50)

di.lax.sub<-droplevels(subset(di.lax, resid(dl3) < cut))
dl4<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = di.lax.sub)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -3.7e+00
qqPlot(resid(dl4))
## [1] 406 438
abline(h = cut <- -50)

di.lax.sub<-droplevels(subset(di.lax.sub, resid(dl4) > cut))
dl5<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position|word), data = di.lax.sub)
## boundary (singular) fit: see ?isSingular
qqPlot(resid(dl5))

## [1] 388 392
plot(allEffects(dl5))

summary(dl5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sylposdur ~ position * sylpos + (1 | subject) + (1 | word) +  
##     (0 + position | word)
##    Data: di.lax.sub
## 
## REML criterion at convergence: 6437.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5433 -0.6713  0.0178  0.6228  3.2631 
## 
## Random effects:
##  Groups   Name                  Variance Std.Dev. Corr 
##  subject  (Intercept)            43.9314  6.628        
##  word     (Intercept)             0.0000  0.000        
##  word.1   positionphrase-medial   1.2255  1.107        
##           positionphrase-final    0.2851  0.534   -1.00
##  Residual                       254.8157 15.963        
## Number of obs: 769, groups:  subject, 10; word, 2
## 
## Fixed effects:
##                               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                    45.0320     2.7581  14.6886  16.327 8.04e-11 ***
## positionphrase-final            8.3210     2.5592   6.3222   3.251 0.016186 *  
## sylposN2                        6.3133     2.2864 751.0147   2.761 0.005898 ** 
## sylposO1                       22.6979     2.2804 751.0022   9.953  < 2e-16 ***
## sylposO2                       35.4649     2.2804 751.0022  15.552  < 2e-16 ***
## positionphrase-final:sylposN2  86.5637     3.2979 751.1636  26.248  < 2e-16 ***
## positionphrase-final:sylposO1   0.2655     3.2250 751.0022   0.082 0.934405    
## positionphrase-final:sylposO2  12.1154     3.2250 751.0022   3.757 0.000185 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- sylpN2 sylpO1 sylpO2 ps-:N2 ps-:O1
## pstnphrs-fn -0.497                                          
## sylposN2    -0.412  0.444                                   
## sylposO1    -0.413  0.446  0.499                            
## sylposO2    -0.413  0.446  0.499  0.500                     
## pstnphr-:N2  0.286 -0.616 -0.693 -0.346 -0.346              
## pstnphr-:O1  0.292 -0.630 -0.353 -0.707 -0.354  0.489       
## pstnphr-:O2  0.292 -0.630 -0.353 -0.354 -0.707  0.489  0.500
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lattice::dotplot(ranef(dl5))
## $subject

## 
## $word

1.4.4 Post-hoc tests (multiple comparisons, Tukey-corrected):

dl5em <- emmeans(dl5, pairwise~position*sylpos)
dl5em.tab<-dl5em$contrasts %>%
  summary(infer = T)
dl5em.tab$p.value<-lapply(dl5em.tab$p.value, FUN = function (x) {
  # levels of p-values
  if (x>=0.05)  p <- "n.\\,s." # not significant
  if (x<0.05)  p <- "<0,05"  # significant
  if (x<=0.01)  p <- "<0,01"  # significant
  if (x<=0.001) p <- "<0,001" # significant
  return(p)
})
dl5em.tab<-as.data.frame(dl5em.tab)
colnames(dl5em.tab)<-c("Comparison","$\\beta$","s","df","Lower CI","Upper CI","t","p")
dl5em.tab.sig<-subset(dl5em.tab, p!="n.\\,s.")
dl5em.tab.sig
##                             Comparison   $\\beta$        s         df
## 3   phrase-medial,N1 - phrase-final,N2 -101.19806 2.644951   6.811796
## 4  phrase-medial,N1 - phrase-medial,O1  -22.69792 2.280421 750.000456
## 5   phrase-medial,N1 - phrase-final,O1  -31.28447 2.559851   5.977351
## 6  phrase-medial,N1 - phrase-medial,O2  -35.46492 2.280421 750.000456
## 7   phrase-medial,N1 - phrase-final,O2  -55.90139 2.559851   5.977351
## 9    phrase-final,N1 - phrase-final,N2  -92.87703 2.377012 750.368370
## 10  phrase-final,N1 - phrase-medial,O1  -14.37688 2.559851   5.977351
## 11   phrase-final,N1 - phrase-final,O1  -22.96344 2.280421 750.000456
## 12  phrase-final,N1 - phrase-medial,O2  -27.14389 2.559851   5.977351
## 13   phrase-final,N1 - phrase-final,O2  -47.58036 2.280421 750.000456
## 14  phrase-medial,N2 - phrase-final,N2  -94.88477 2.649853   6.863926
## 15 phrase-medial,N2 - phrase-medial,O1  -16.38463 2.286408 750.016868
## 16  phrase-medial,N2 - phrase-final,O1  -24.97118 2.564669   6.024249
## 17 phrase-medial,N2 - phrase-medial,O2  -29.15164 2.286408 750.016868
## 18  phrase-medial,N2 - phrase-final,O2  -49.58810 2.564669   6.024249
## 19  phrase-final,N2 - phrase-medial,O1   78.50015 2.644951   6.811796
## 20   phrase-final,N2 - phrase-final,O1   69.91359 2.377012 750.368370
## 21  phrase-final,N2 - phrase-medial,O2   65.73314 2.644951   6.811796
## 22   phrase-final,N2 - phrase-final,O2   45.29667 2.377012 750.368370
## 24 phrase-medial,O1 - phrase-medial,O2  -12.76701 2.280421 750.000456
## 25  phrase-medial,O1 - phrase-final,O2  -33.20348 2.559851   5.977351
## 27   phrase-final,O1 - phrase-final,O2  -24.61692 2.280421 750.000456
## 28  phrase-medial,O2 - phrase-final,O2  -20.43647 2.559851   5.977351
##      Lower CI   Upper CI          t      p
## 3  -112.16716 -90.228969 -38.260840 <0,001
## 4   -29.62924 -15.766594  -9.953389 <0,001
## 5   -42.38136 -20.187581 -12.221210 <0,001
## 6   -42.39625 -28.533601 -15.551920 <0,001
## 7   -66.99828 -44.804503 -21.837756 <0,001
## 9  -100.10193 -85.652128 -39.073021 <0,001
## 10  -25.47377  -3.279992  -5.616297  <0,05
## 11  -29.89476 -16.032113 -10.069823 <0,001
## 12  -38.24078 -16.046999 -10.603700 <0,001
## 13  -54.51168 -40.649035 -20.864726 <0,001
## 14 -105.84781 -83.921738 -35.807557 <0,001
## 15  -23.33415  -9.435109  -7.166100 <0,001
## 16  -36.05765 -13.884715  -9.736609 <0,001
## 17  -36.10116 -22.202115 -12.749970 <0,001
## 18  -60.67457 -38.501636 -19.335086 <0,001
## 19   67.53105  89.469239  29.679240 <0,001
## 20   62.68869  77.138491  29.412389 <0,001
## 21   54.76405  76.702232  24.852305 <0,001
## 22   38.07177  52.521570  19.056141 <0,001
## 24  -19.69833  -5.835684  -5.598531 <0,001
## 25  -44.30036 -22.106586 -12.970865 <0,001
## 27  -31.54824 -17.685599 -10.794903 <0,001
## 28  -31.53336  -9.339580  -7.983462  <0,01

1.4.5 Disyllabic effect plot

Preparation

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

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

1.5 Figure 3: Mono- and disyllabic effect plot

ggpubr::ggarrange(ef1b.plot, ef2b.plot, common.legend = T,
                  labels = c("a)","b)"),font.label = list(size=11,face="plain"),align = "h",nrow=1)

ggsave(file="mono-di-acoustic-effects.pdf", width = 170, height = 100, units = "mm")

1.6 Table 3 and 4

x<-capture.output(
   texreg::texreg(list(texreg::extract(m5, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                       # texreg::extract(d5, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                       # texreg::extract(mf3.1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
   ),
   digits = 1,booktabs = T,custom.model.names = c("Duration in monosyllabic words"),
   float.pos = "!htbp", label = "tab:acousticresults.mono", use.packages = F,single.row = T,caption.above = T,
   caption="Model predictions for the acoustic duration of syllable position in  monosyllabic words
   (standards error in parentheses) for the fixed effects  phrasal position (phrase-medial vs. phrase-final),
   syllable constituent (\\emph{sylcon}, O = onset, N = nucleus, C = coda), and tenseness (tense vs. lax).",
   dcolumn = T, fontsize = "scriptsize",center=F)
 )
 x[5]<-paste0(x[5], "\\centering")
 x[27]<-paste0(x[27],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(m5.r2dummy)[1],2),"/",round(r.squaredGLMM(m5.r2dummy)[2],2),
               # "&", round(r.squaredGLMM(mf3.1)[1],2), "/", round(r.squaredGLMM(mf3.1)[2],2),
               "\\\\ \n")
 cat(x, sep = "\n")
## 
## \begin{table}[!htbp]
## \caption{Model predictions for the acoustic duration of syllable position in  monosyllabic words
##    (standards error in parentheses) for the fixed effects  phrasal position (phrase-medial vs. phrase-final),
##    syllable constituent (\emph{sylcon}, O = onset, N = nucleus, C = coda), and tenseness (tense vs. lax).}\centering
## \begin{scriptsize}
## \begin{tabular}{l D{)}{)}{10)3} }
## \toprule
##  & \multicolumn{1}{c}{Duration in monosyllabic words} \\
## \midrule
## (Intercept)                               & 93.8 \; (6.3)^{***}  \\
## positionphrase-final                      & 10.6 \; (4.7)^{*}    \\
## sylposN                                   & 6.3 \; (3.3)         \\
## sylposC                                   & -29.4 \; (3.3)^{***} \\
## tensenesslax                              & -5.9 \; (8.1)        \\
## positionphrase-final:sylposN              & 34.1 \; (4.6)^{***}  \\
## positionphrase-final:sylposC              & 53.2 \; (4.6)^{***}  \\
## positionphrase-final:tensenesslax         & -2.7 \; (5.4)        \\
## sylposN:tensenesslax                      & -31.5 \; (4.6)^{***} \\
## sylposC:tensenesslax                      & 13.2 \; (4.6)^{**}   \\
## positionphrase-final:sylposN:tensenesslax & -16.1 \; (6.5)^{*}   \\
## positionphrase-final:sylposC:tensenesslax & 11.9 \; (6.5)        \\
## \midrule
## AIC                                       & 22939.7              \\
## Num. obs.                                 & 2342                 \\
## Num. groups: subject                      & 10                   \\
## Num. groups: word                         & 8                    \\\midrule R$_m^2$/R$_c^2$ & 0.29/0.41\\ 
## 
## \bottomrule
## \multicolumn{2}{l}{\tiny{$^{***}p<0.001$, $^{**}p<0.01$, $^*p<0.05$}}
## \end{tabular}
## \end{scriptsize}
## \label{tab:acousticresults.mono}
## \end{table}
 x<-capture.output(
   texreg::texreg(list(texreg::extract(d5, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                       # texreg::extract(d5, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                       # texreg::extract(mf3.1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
   ),
   digits = 1,booktabs = T,custom.model.names = c("Duration in disyllabic words"),
   float.pos = "!htbp", label = "tab:acousticresults.di", use.packages = F,single.row = T,caption.above = T,
   caption="Model predictions for the acoustic duration of syllable position in disyllabic words
   (standards error in parentheses) for the fixed effects  phrasal position (phrase-medial vs. phrase-final),
   syllable position (\\emph{sylcon}, O1/O2 = first/second onset, N1/N2 = first/second nucleus), and tenseness (tense vs. lax).",
   dcolumn = T, fontsize = "scriptsize",center=F)
 )
 x[5]<-paste0(x[5], "\\centering")
 x[31]<-paste0(x[31],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(d5.r2dummy)[1],2),"/", round(r.squaredGLMM(d5.r2dummy)[2],2),
               "\\\\ \n")
 cat(x, sep = "\n")
## 
## \begin{table}[!htbp]
## \caption{Model predictions for the acoustic duration of syllable position in disyllabic words
##    (standards error in parentheses) for the fixed effects  phrasal position (phrase-medial vs. phrase-final),
##    syllable position (\emph{sylcon}, O1/O2 = first/second onset, N1/N2 = first/second nucleus), and tenseness (tense vs. lax).}\centering
## \begin{scriptsize}
## \begin{tabular}{l D{)}{)}{11)3} }
## \toprule
##  & \multicolumn{1}{c}{Duration in disyllabic words} \\
## \midrule
## (Intercept)                                & 75.7 \; (14.4)^{***} \\
## positionphrase-final                       & 4.7 \; (2.9)         \\
## sylposN1                                   & -1.7 \; (18.9)       \\
## sylposO2                                   & 3.8 \; (12.5)        \\
## sylposN2                                   & -31.9 \; (18.4)      \\
## tensenesslax                               & -7.8 \; (19.9)       \\
## positionphrase-final:sylposN1              & 14.5 \; (3.1)^{***}  \\
## positionphrase-final:sylposO2              & 11.8 \; (3.1)^{***}  \\
## positionphrase-final:sylposN2              & 100.0 \; (3.1)^{***} \\
## positionphrase-final:tensenesslax          & 3.4 \; (4.1)         \\
## sylposN1:tensenesslax                      & -21.3 \; (26.8)      \\
## sylposO2:tensenesslax                      & 8.8 \; (17.7)        \\
## sylposN2:tensenesslax                      & 15.3 \; (26.1)       \\
## positionphrase-final:sylposN1:tensenesslax & -14.1 \; (4.4)^{**}  \\
## positionphrase-final:sylposO2:tensenesslax & 0.4 \; (4.4)         \\
## positionphrase-final:sylposN2:tensenesslax & -11.1 \; (4.5)^{*}   \\
## \midrule
## AIC                                        & 12857.9              \\
## Num. obs.                                  & 1541                 \\
## Num. groups: subject                       & 10                   \\
## Num. groups: word                          & 4                    \\\midrule R$_m^2$/R$_c^2$ & 0.68/0.81\\ 
## 
## \bottomrule
## \multicolumn{2}{l}{\tiny{$^{***}p<0.001$, $^{**}p<0.01$, $^*p<0.05$}}
## \end{tabular}
## \end{scriptsize}
## \label{tab:acousticresults.di}
## \end{table}

1.7 Table 9 and 10: Post-hoc tests (multiple comparisons, Tukey-corrected):

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

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

print(xtable::xtable(m5em.tab.sig,label="tab:acoustic.mono.post",
                      caption = "Significant contrasts for monosyllabic target words with duration as dependent variable,
                      the predictors position (phrase-medial/-final), syllable position (O = Onset, N = Nucleus, C = Coda),
                      and tenseness (tense/lax), as well as  estimated value $\\beta$, standard deviation (s), degrees of freedom (df), 95\\,\\% confidence
                      intervals, and t-values, Tukey-corrected for  multiple comparisons.",
                      align = "llrrrrrrr",digits = c(0,0,1,1,1,2,2,2,2)),
       include.colnames = T,include.rownames=T,booktabs = T, size = "\\scriptsize",
       caption.placement = "top", sanitize.text.function = function(x) x)
## % latex table generated in R 3.6.0 by xtable 1.8-4 package
## % Wed Jan 20 16:03:46 2021
## \begin{table}[ht]
## \centering
## \caption{Significant contrasts for monosyllabic target words with duration as dependent variable,
##                       the predictors position (phrase-medial/-final), syllable position (O = Onset, N = Nucleus, C = Coda),
##                       and tenseness (tense/lax), as well as  estimated value $\beta$, standard deviation (s), degrees of freedom (df), 95\,\% confidence
##                       intervals, and t-values, Tukey-corrected for  multiple comparisons.} 
## \label{tab:acoustic.mono.post}
## \begingroup\scriptsize
## \begin{tabular}{llrrrrrrr}
##   \toprule
##  & Comparison & $\beta$ & s & df & Lower CI & Upper CI & t & p \\ 
##   \midrule
## 1 & phrase-medial,N,tense - phrase-final,N,tense & -44.7 & 4.6 & 25.9 & -61.31 & -28.10 & -9.64 & <0.001 \\ 
##   2 & phrase-medial,N,tense - phrase-medial,C,tense & 35.7 & 3.3 & 2291.0 & 25.00 & 46.38 & 10.92 & <0.001 \\ 
##   3 & phrase-medial,N,tense - phrase-final,C,tense & -28.0 & 4.6 & 25.9 & -44.63 & -11.42 & -6.04 & <0.001 \\ 
##   4 & phrase-medial,N,tense - phrase-medial,N,lax & 37.4 & 8.1 & 7.6 & 1.54 & 73.24 & 4.63 & <0.05 \\ 
##   5 & phrase-medial,N,tense - phrase-final,C,lax & -44.5 & 9.7 & 8.1 & -86.77 & -2.14 & -4.57 & <0.05 \\ 
##   6 & phrase-final,N,tense - phrase-medial,C,tense & 80.4 & 4.6 & 25.9 & 63.79 & 97.00 & 17.34 & <0.001 \\ 
##   7 & phrase-final,N,tense - phrase-final,C,tense & 16.7 & 3.2 & 2291.0 & 6.15 & 27.21 & 5.18 & <0.001 \\ 
##   8 & phrase-final,N,tense - phrase-medial,N,lax & 82.1 & 10.0 & 8.9 & 39.74 & 124.45 & 8.23 & <0.001 \\ 
##   9 & phrase-final,N,tense - phrase-final,N,lax & 56.2 & 10.7 & 6.8 & 7.03 & 105.32 & 5.24 & <0.05 \\ 
##   10 & phrase-final,N,tense - phrase-medial,C,lax & 73.1 & 10.0 & 8.9 & 30.75 & 115.46 & 7.33 & <0.01 \\ 
##   11 & phrase-medial,C,tense - phrase-final,C,tense & -63.7 & 4.6 & 25.9 & -80.32 & -47.11 & -13.74 & <0.001 \\ 
##   12 & phrase-medial,C,tense - phrase-final,C,lax & -80.1 & 9.7 & 8.1 & -122.47 & -37.83 & -8.24 & <0.001 \\ 
##   13 & phrase-final,C,tense - phrase-medial,N,lax & 65.4 & 10.0 & 8.9 & 23.06 & 107.77 & 6.56 & <0.01 \\ 
##   14 & phrase-final,C,tense - phrase-medial,C,lax & 56.4 & 10.0 & 8.9 & 14.07 & 98.78 & 5.66 & <0.01 \\ 
##   15 & phrase-medial,N,lax - phrase-final,N,lax & -25.9 & 4.6 & 25.6 & -42.50 & -9.34 & -5.60 & <0.001 \\ 
##   16 & phrase-medial,N,lax - phrase-final,C,lax & -81.8 & 4.6 & 25.6 & -98.42 & -65.27 & -17.70 & <0.001 \\ 
##   17 & phrase-final,N,lax - phrase-medial,C,lax & 16.9 & 4.6 & 25.7 & 0.34 & 33.51 & 3.66 & <0.05 \\ 
##   18 & phrase-final,N,lax - phrase-final,C,lax & -55.9 & 3.2 & 2291.0 & -66.43 & -45.43 & -17.42 & <0.001 \\ 
##   19 & phrase-medial,C,lax - phrase-final,C,lax & -72.9 & 4.6 & 25.7 & -89.44 & -56.27 & -15.74 & <0.001 \\ 
##    \bottomrule
## \end{tabular}
## \endgroup
## \end{table}
 print(xtable::xtable(d5em.tab.sig,label="tab:acoustic.di.post",
                      caption = "Significant contrasts for disyllabic target words with duration as dependent variable,
                      the predictors position (phrase-medial/-final), syllable position (O1/O2 = first/second onset, N1/N2 = first/second nucleus, C = coda),
                      and tenseness (tense/lax), as well as  estimated value $\\beta$, standard deviation (s), degrees of freedom (df), 95\\,\\% confidence
                      intervals, and t-values, Tukey-corrected for  multiple comparisons.",
                      align = "llrrrrrrr",digits = c(0,0,1,1,1,2,2,2,2)),
       include.colnames = T,include.rownames=T,booktabs = T, size = "\\scriptsize",
       caption.placement = "top", sanitize.text.function = function(x) x)
## % latex table generated in R 3.6.0 by xtable 1.8-4 package
## % Wed Jan 20 16:03:46 2021
## \begin{table}[ht]
## \centering
## \caption{Significant contrasts for disyllabic target words with duration as dependent variable,
##                       the predictors position (phrase-medial/-final), syllable position (O1/O2 = first/second onset, N1/N2 = first/second nucleus, C = coda),
##                       and tenseness (tense/lax), as well as  estimated value $\beta$, standard deviation (s), degrees of freedom (df), 95\,\% confidence
##                       intervals, and t-values, Tukey-corrected for  multiple comparisons.} 
## \label{tab:acoustic.di.post}
## \begingroup\scriptsize
## \begin{tabular}{llrrrrrrr}
##   \toprule
##  & Comparison & $\beta$ & s & df & Lower CI & Upper CI & t & p \\ 
##   \midrule
## 1 & phrase-medial,N1,tense - phrase-final,N1,tense & -19.2 & 2.9 & 6.0 & -34.02 & -4.32 & -6.60 & <0.05 \\ 
##   2 & phrase-medial,N1,tense - phrase-medial,N2,tense & 30.3 & 2.4 & 6.2 & 18.36 & 42.17 & 12.88 & <0.001 \\ 
##   3 & phrase-medial,N1,tense - phrase-final,N2,tense & -74.4 & 2.7 & 2.9 & -95.43 & -53.38 & -27.51 & <0.01 \\ 
##   4 & phrase-final,N1,tense - phrase-medial,N2,tense & 49.4 & 3.3 & 2.5 & 20.01 & 78.86 & 14.82 & <0.05 \\ 
##   5 & phrase-final,N1,tense - phrase-final,N2,tense & -55.2 & 2.4 & 6.4 & -67.13 & -43.34 & -23.24 & <0.001 \\ 
##   6 & phrase-medial,O2,tense - phrase-final,O2,tense & -16.5 & 2.9 & 6.0 & -31.36 & -1.66 & -5.69 & <0.05 \\ 
##   7 & phrase-medial,O2,tense - phrase-final,O2,lax & -21.4 & 3.0 & 5.6 & -37.12 & -5.70 & -7.20 & <0.05 \\ 
##   8 & phrase-medial,O2,tense - phrase-final,N2,lax & -68.8 & 7.1 & 2.3 & -135.86 & -1.71 & -9.70 & <0.05 \\ 
##   9 & phrase-final,O2,tense - phrase-medial,N2,tense & 52.2 & 4.6 & 2.2 & 6.55 & 97.90 & 11.29 & <0.05 \\ 
##   10 & phrase-final,O2,tense - phrase-medial,N1,lax & 51.1 & 5.5 & 2.6 & 3.55 & 98.55 & 9.23 & <0.05 \\ 
##   11 & phrase-final,O2,tense - phrase-medial,O2,lax & 15.5 & 2.8 & 6.2 & 1.44 & 29.51 & 5.59 & <0.05 \\ 
##   12 & phrase-final,O2,tense - phrase-medial,N2,lax & 44.7 & 5.1 & 2.7 & 2.60 & 86.80 & 8.76 & <0.05 \\ 
##   13 & phrase-medial,N2,tense - phrase-final,N2,tense & -104.7 & 2.9 & 6.4 & -119.44 & -89.90 & -35.50 & <0.001 \\ 
##   14 & phrase-medial,N2,tense - phrase-final,O2,lax & -57.1 & 4.9 & 2.8 & -95.74 & -18.50 & -11.74 & <0.05 \\ 
##   15 & phrase-medial,N2,tense - phrase-final,N2,lax & -104.5 & 8.1 & 2.2 & -184.38 & -24.62 & -12.96 & <0.05 \\ 
##   16 & phrase-final,N2,tense - phrase-medial,N1,lax & 103.5 & 8.3 & 2.3 & 23.98 & 183.02 & 12.53 & <0.05 \\ 
##   17 & phrase-final,N2,tense - phrase-medial,O2,lax & 67.9 & 6.7 & 2.4 & 4.74 & 131.11 & 10.09 & <0.05 \\ 
##   18 & phrase-final,N2,tense - phrase-medial,N2,lax & 97.1 & 8.0 & 2.2 & 17.49 & 176.81 & 12.17 & <0.05 \\ 
##   19 & phrase-medial,N1,lax - phrase-final,O2,lax & -56.0 & 5.1 & 2.2 & -106.93 & -4.97 & -11.07 & <0.05 \\ 
##   20 & phrase-medial,N1,lax - phrase-final,N2,lax & -103.3 & 2.7 & 3.1 & -123.62 & -83.03 & -37.66 & <0.001 \\ 
##   21 & phrase-final,N1,lax - phrase-final,N2,lax & -94.8 & 2.4 & 7.0 & -106.60 & -83.09 & -39.02 & <0.001 \\ 
##   22 & phrase-medial,O2,lax - phrase-final,O2,lax & -20.4 & 2.9 & 6.0 & -35.24 & -5.52 & -7.02 & <0.05 \\ 
##   23 & phrase-final,O2,lax - phrase-medial,N2,lax & 49.6 & 4.6 & 2.2 & 3.91 & 95.29 & 10.72 & <0.05 \\ 
##   24 & phrase-medial,N2,lax - phrase-final,N2,lax & -97.0 & 3.0 & 6.7 & -111.64 & -82.31 & -32.51 & <0.001 \\ 
##    \bottomrule
## \end{tabular}
## \endgroup
## \end{table}

2 Articulatory analysis

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

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


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

2.1 Data removal and data used

Remove all items with closing or opening displacement below 1mm Monosyllabic 389 target words were originally included:

 addmargins(table(mono.art$word,mono.art$position))
##        
##         phrase-medial phrase-final Sum
##   Bahn             50           49  99
##   Bann             48           51  99
##   Beet             47           51  98
##   Bett             48           49  97
##   Ruhm             49           47  96
##   Rum              49           48  97
##   Stiel            48           51  99
##   still            50           49  99
##   Sum             389          395 784

Disyllabic 198 target words were originally included:

 addmargins(table(di.art$word,di.art$position))
##         
##          phrase-medial phrase-final Sum
##   Huete             49           50  99
##   Huette            51           48  99
##   Miete             49           47  96
##   Mitte             49           50  99
##   Sum              198          195 393
 ggplot(mono.art, aes(x=word, y=closdisp))+
   geom_point()+
   geom_hline(yintercept=1)

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

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

 mono.art.norm<- mono.art[mono.art$opdisp>1 & mono.art$closdisp>1,]
 nrow(mono.art)-nrow(mono.art.norm)
## [1] 36
 #' 748 stimuli after data exclusion:
 addmargins(table(mono.art.norm$word,mono.art.norm$position))
##        
##         phrase-medial phrase-final Sum
##   Bahn             50           49  99
##   Bann             48           51  99
##   Beet             47           51  98
##   Bett             48           48  96
##   Ruhm             39           33  72
##   Rum              49           40  89
##   Stiel            47           51  98
##   still            49           48  97
##   Sum             377          371 748
 # nrow(di.art[di.art$opdisp<1 | di.art$closdisp<1 ,])
 di.art.norm<- di.art[di.art$opdisp>1 & di.art$closdisp>1,]
 nrow(di.art)-nrow(di.art.norm)
## [1] 7
 # invalid data of closing duration outlier, remove:
 di.art.norm<-di.art.norm[di.art.norm$closdur>15,]
 nrow(di.art)-nrow(di.art.norm)
## [1] 8

385 stimuli after data exclusion:

 addmargins(table(di.art.norm$word,di.art.norm$position))
##         
##          phrase-medial phrase-final Sum
##   Huete             49           49  98
##   Huette            49           45  94
##   Miete             49           46  95
##   Mitte             48           50  98
##   Sum              195          190 385

Total data excluded in the articulatory analysis

 # nrow(mono.art)/nrow(mono.art.norm)-1

Monosyllabic

 sprintf("%1.2f%%",round((1- nrow(mono.art.norm)/800)*100),2)
## [1] "6.00%"

Disyllabic

 #nrow(di.art)/nrow(di.art.norm)-1
 sprintf("%1.2f%%",round((1- nrow(di.art.norm)/400)*100),2)
## [1] "4.00%"

** Total total **

 sprintf("%1.2f%%",round((1- ((nrow(di.art.norm)+nrow(mono.art.norm))/(400+800)) )*100,2))
## [1] "5.58%"

Number of target words used in the articulatory analysis

 addmargins(table(rbind(mono.art.norm,di.art.norm)$word, rbind(mono.art.norm,di.art.norm)$position))
##         
##          phrase-medial phrase-final  Sum
##   Bahn              50           49   99
##   Bann              48           51   99
##   Beet              47           51   98
##   Bett              48           48   96
##   Ruhm              39           33   72
##   Rum               49           40   89
##   Stiel             47           51   98
##   still             49           48   97
##   Huete             49           49   98
##   Huette            49           45   94
##   Miete             49           46   95
##   Mitte             48           50   98
##   Sum              572          561 1133
 mono.art$wordpair<-NULL
 mono.art$wordpair[mono.art$word=="Bahn" | mono.art$word=="Bann"] <-"Bahn/Bann"
 mono.art$wordpair[mono.art$word=="Beet" | mono.art$word=="Bett"] <-"Beet/Bett"
 mono.art$wordpair[mono.art$word=="Ruhm" | mono.art$word=="Rum"] <-"Ruhm/Rum"
 mono.art$wordpair[mono.art$word=="Stiel" | mono.art$word=="still"] <-"Stiel/still"


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

2.2 Monosyllabic target words

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

2.2.1 Closing duration

 # Closing duration mono ####
 clodur1<-lmer(log(closdur) ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art)
 qqPlot(resid(clodur1))
## 11312  2532 
##   674   740
 abline(h = cut <- .6)

 mono.art.sub1<-droplevels(subset(mono.art, resid(clodur1) < cut))
 clodur2<-lmer(log(closdur) ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art.sub1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00228399 (tol = 0.002, component 1)
 qqPlot(resid(clodur2))

## 217 220 
##  86  89
 summary(clodur2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(closdur) ~ position + tenseness + (1 + position + tenseness |  
##     subject) + (1 | word)
##    Data: mono.art.sub1
## 
## REML criterion at convergence: 241.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7846 -0.5867  0.0197  0.6400  2.6637 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr       
##  subject  (Intercept)          0.008854 0.09410             
##           positionphrase-final 0.002446 0.04945  -0.20      
##           tensenesslax         0.008308 0.09115  -0.13  0.23
##  word     (Intercept)          0.029372 0.17138             
##  Residual                      0.072934 0.27006             
## Number of obs: 760, groups:  subject, 10; word, 8
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)           3.84979    0.09230  7.56129  41.710 3.25e-10 ***
## positionphrase-final  0.32404    0.02509  9.10505  12.915 3.67e-07 ***
## tensenesslax         -0.06291    0.12610  6.66174  -0.499    0.634    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.123       
## tensenesslx -0.657  0.034
## convergence code: 0
## Model failed to converge with max|grad| = 0.00228399 (tol = 0.002, component 1)
 plot(allEffects(clodur2))

 clodur3.int<-lmer(log(closdur) ~ position  * tenseness + (1|subject) + (1|word), data = mono.art.sub1)
 anova(clodur2, clodur3.int) # not better
## refitting model(s) with ML (instead of REML)
## Data: mono.art.sub1
## Models:
## clodur3.int: log(closdur) ~ position * tenseness + (1 | subject) + (1 | word)
## clodur2: log(closdur) ~ position + tenseness + (1 + position + tenseness | 
## clodur2:     subject) + (1 | word)
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## clodur3.int  7 254.00 286.44 -120.00   240.00                           
## clodur2     11 251.44 302.41 -114.72   229.44 10.565      4    0.03192 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 # failed

 clodur2.ms<-lmer(closdur ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art.sub1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0351818 (tol = 0.002, component 1)
 summary(clodur2.ms)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdur ~ position + tenseness + (1 + position + tenseness |  
##     subject) + (1 | word)
##    Data: mono.art.sub1
## 
## REML criterion at convergence: 6399.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5627 -0.6386 -0.0711  0.4827  3.9194 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr       
##  subject  (Intercept)           22.28    4.720              
##           positionphrase-final  10.23    3.199    0.84      
##           tensenesslax          19.35    4.399   -0.15  0.23
##  word     (Intercept)           83.25    9.124              
##  Residual                      250.50   15.827              
## Number of obs: 760, groups:  subject, 10; word, 8
## 
## Fixed effects:
##                      Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)            49.875      4.903  7.571  10.173 1.11e-05 ***
## positionphrase-final   18.331      1.531  9.470  11.970 4.96e-07 ***
## tensenesslax           -3.672      6.699  6.671  -0.548    0.601    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn  0.081       
## tensenesslx -0.663  0.032
## convergence code: 0
## Model failed to converge with max|grad| = 0.0351818 (tol = 0.002, component 1)
 plot(allEffects(clodur2.ms))

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

## 12910  1211 
##   258   200
 plot(allEffects(cd.lax.mono.1))

 summary(cd.lax.mono.1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(closdur) ~ position + (1 | subject) + (1 | word)
##    Data: lax.cd.mono
## 
## REML criterion at convergence: 98.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2177 -0.6486 -0.0523  0.6784  2.7575 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 0.01569  0.1252  
##  word     (Intercept) 0.02665  0.1633  
##  Residual             0.06801  0.2608  
## Number of obs: 380, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                       Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            3.78426    0.09267   4.62551   40.84 4.27e-07 ***
## positionphrase-final   0.32868    0.02680 366.23567   12.27  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pstnphrs-fn -0.143

2.2.2 Plateau duration

 platdur1<-lmer(log(platdur) ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art)
## boundary (singular) fit: see ?isSingular
 qqPlot(resid(platdur1))
## 745 546 
## 259 188
 abline(h = cut <- c(-1.1,1.3))

 mono.art.plat1<-droplevels(subset(mono.art, resid(platdur1) > cut[1] & resid(platdur1) < cut[2]))
 platdur2<-lmer(log(platdur) ~ position  + tenseness + (1|subject) + (1|word), data = mono.art.plat1)
 qqPlot(resid(platdur2))

##  219 1018 
##   88  331
 summary(platdur2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(platdur) ~ position + tenseness + (1 | subject) + (1 | word)
##    Data: mono.art.plat1
## 
## REML criterion at convergence: 1040.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2799 -0.6406  0.0006  0.6190  3.0363 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 0.02695  0.1642  
##  word     (Intercept) 0.10009  0.3164  
##  Residual             0.21825  0.4672  
## Number of obs: 748, groups:  subject, 10; word, 8
## 
## Fixed effects:
##                       Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            3.43041    0.16905   7.40951  20.292 9.17e-08 ***
## positionphrase-final   0.81460    0.03426 730.53104  23.779  < 2e-16 ***
## tensenesslax           0.20352    0.22630   6.00295   0.899    0.403    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.096       
## tensenesslx -0.669 -0.001
 plot(allEffects(platdur2))

 platdur3<-lmer(log(platdur) ~ position  * tenseness + (1|subject) + (1|word), data = mono.art.plat1)
 anova(platdur2,platdur3) #no
## refitting model(s) with ML (instead of REML)
## Data: mono.art.plat1
## Models:
## platdur2: log(platdur) ~ position + tenseness + (1 | subject) + (1 | word)
## platdur3: log(platdur) ~ position * tenseness + (1 | subject) + (1 | word)
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## platdur2  6 1044.3 1072.0 -516.14   1032.3                         
## platdur3  7 1045.6 1077.9 -515.81   1031.6 0.6731      1      0.412
 platdur2.ms<-lmer(platdur ~ position  + tenseness + (1|subject) + (1|word), data = mono.art.plat1)
 summary(platdur2.ms)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: platdur ~ position + tenseness + (1 | subject) + (1 | word)
##    Data: mono.art.plat1
## 
## REML criterion at convergence: 7864.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9226 -0.6470 -0.1173  0.3996  5.7967 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept)  207.5   14.41   
##  word     (Intercept)  747.3   27.34   
##  Residual             2081.8   45.63   
## Number of obs: 748, groups:  subject, 10; word, 8
## 
## Fixed effects:
##                      Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)            29.878     14.690   7.480   2.034   0.0789 .  
## positionphrase-final   62.360      3.346 730.642  18.640   <2e-16 ***
## tensenesslax           14.227     19.617   6.007   0.725   0.4956    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.108       
## tensenesslx -0.668 -0.002

2.2.3 Closing displacement

 closdis1<-lmer(closdisp ~ position  + tenseness +(1+position+tenseness|subject) + (1|word), data = mono.art)
## boundary (singular) fit: see ?isSingular
 qqPlot(resid(closdis1))
## 1384 1792 
##  467  690
 abline(h = cut <- 10)

 mono.art.disp1<-droplevels(subset(mono.art, resid(closdis1) < cut))
 closdis2<-lmer(closdisp ~ position  + tenseness +(1+position|subject) + (1|word), data = mono.art.disp1)
 qqPlot(resid(closdis2))

## 1792  378 
##  689  135
 summary(closdis2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdisp ~ position + tenseness + (1 + position | subject) +  
##     (1 | word)
##    Data: mono.art.disp1
## 
## REML criterion at convergence: 3526.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7968 -0.6069 -0.0275  0.5838  3.3944 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr 
##  subject  (Intercept)           2.0615  1.4358        
##           positionphrase-final  0.3725  0.6104   -0.72
##  word     (Intercept)          32.7455  5.7224        
##  Residual                       4.7875  2.1880        
## Number of obs: 783, groups:  subject, 10; word, 8
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)           7.64946    2.90017  6.30870   2.638   0.0369 *
## positionphrase-final  0.62736    0.24847  8.91719   2.525   0.0327 *
## tensenesslax         -0.05315    4.04935  5.99701  -0.013   0.9900  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.104       
## tensenesslx -0.698  0.000
 plot(allEffects(closdis2))

 # final is more displaced than medial

2.2.4 Closing velocity

 closvel1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art)
 qqPlot(resid(closvel1))
## 1024 1030 
##  345  351
 abline(h = cut <- 18)

 mono.art.vel1<-droplevels(subset(mono.art, resid(closvel1) < cut))
 closvel2<-lmer(PVEL..V.cm.sec. ~ position  + tenseness + (1+position|subject) + (1|word), data = mono.art.vel1)
 qqPlot(resid(closvel2))

## 1792 1384 
##  688  465
 summary(closvel2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + tenseness + (1 + position | subject) +  
##     (1 | word)
##    Data: mono.art.vel1
## 
## REML criterion at convergence: 4920.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3003 -0.6261  0.0040  0.5993  3.2141 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr 
##  subject  (Intercept)           27.744   5.267        
##           positionphrase-final   5.859   2.421   -0.61
##  word     (Intercept)          162.302  12.740        
##  Residual                       28.321   5.322        
## Number of obs: 782, groups:  subject, 10; word, 8
## 
## Fixed effects:
##                      Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)            20.466      6.592  6.837   3.105  0.01774 * 
## positionphrase-final   -3.248      0.855  8.922  -3.799  0.00429 **
## tensenesslax            1.661      9.016  6.001   0.184  0.85995   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.151       
## tensenesslx -0.684  0.000
 plot(allEffects(closvel2))

 # dotplot(ranef(closvel2))
 # final is sign. slower than medial

In none of the models with position and tenseness as fixed and subjects and words as random effects does tenseness affect the outcome variable ## Articulatory disyllabic analysis #### ### Closing duration

 di.clodur1<-lmer(closdur ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = di.art)
## boundary (singular) fit: see ?isSingular
 qqPlot(resid(di.clodur1))
## 323 443 
##  78  98
 abline(h = cut <- 20)

 di.art.sub1<-droplevels(subset(di.art, resid(di.clodur1) < cut))
 di.clodur2<-lmer(closdur ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = di.art.sub1)
 di.clodur3<-lmer(closdur ~ position  + tenseness +
                    (1+position+tenseness|subject) + (1|word), data = di.art.sub1)
 anova(di.clodur2,di.clodur3,refit=F)
## Data: di.art.sub1
## Models:
## di.clodur2: closdur ~ position + tenseness + (1 + position + tenseness | 
## di.clodur2:     subject) + (1 | word)
## di.clodur3: closdur ~ position + tenseness + (1 + position + tenseness | 
## di.clodur3:     subject) + (1 | word)
##            Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## di.clodur2 11 2754.3 2797.4 -1366.1   2732.3                        
## di.clodur3 11 2754.3 2797.4 -1366.1   2732.3     0      0          1
 qqPlot(resid(di.clodur2))

## 625 999 
## 114 217
 summary(di.clodur2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdur ~ position + tenseness + (1 + position + tenseness |  
##     subject) + (1 | word)
##    Data: di.art.sub1
## 
## REML criterion at convergence: 2732.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6575 -0.5757 -0.0672  0.6350  2.5307 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr       
##  subject  (Intercept)          42.93    6.552               
##           positionphrase-final 32.31    5.684     0.51      
##           tensenesslax         53.15    7.290    -0.64 -0.89
##  word     (Intercept)          13.30    3.647               
##  Residual                      78.35    8.852               
## Number of obs: 373, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                      Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)            47.447      3.403  4.856  13.944  4.2e-05 ***
## positionphrase-final    3.382      2.020  8.701   1.674    0.130    
## tensenesslax            2.480      4.412  3.644   0.562    0.607    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn  0.219       
## tensenesslx -0.676 -0.417
 plot(allEffects(di.clodur2))

 di.clodur.int<-lmer(closdur ~ position  * tenseness + (1+position+tenseness|subject) + (1|word), data = di.art.sub1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00514629 (tol = 0.002, component 1)
 anova(di.clodur2,di.clodur.int) # not better
## refitting model(s) with ML (instead of REML)
## Data: di.art.sub1
## Models:
## di.clodur2: closdur ~ position + tenseness + (1 + position + tenseness | 
## di.clodur2:     subject) + (1 | word)
## di.clodur.int: closdur ~ position * tenseness + (1 + position + tenseness | 
## di.clodur.int:     subject) + (1 | word)
##               Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## di.clodur2    11 2765.3 2808.5 -1371.7   2743.3                        
## di.clodur.int 12 2767.3 2814.4 -1371.7   2743.3 1e-04      1     0.9915
 summary(di.clodur.int)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdur ~ position * tenseness + (1 + position + tenseness |  
##     subject) + (1 | word)
##    Data: di.art.sub1
## 
## REML criterion at convergence: 2729.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6529 -0.5757 -0.0678  0.6338  2.5283 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr       
##  subject  (Intercept)          42.88    6.548               
##           positionphrase-final 32.33    5.686     0.51      
##           tensenesslax         53.14    7.289    -0.64 -0.89
##  word     (Intercept)          13.37    3.656               
##  Residual                      78.58    8.865               
## Number of obs: 373, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                                    Estimate Std. Error        df t value
## (Intercept)                        47.44274    3.43796   4.98555  13.800
## positionphrase-final                3.39180    2.23267  12.91091   1.519
## tensenesslax                        2.48910    4.51153   3.91963   0.552
## positionphrase-final:tensenesslax  -0.01866    1.84432 340.87161  -0.010
##                                   Pr(>|t|)    
## (Intercept)                       3.66e-05 ***
## positionphrase-final                 0.153    
## tensenesslax                         0.611    
## positionphrase-final:tensenesslax    0.992    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- tnsnss
## pstnphrs-fn  0.139              
## tensenesslx -0.683 -0.284       
## pstnphrs-f:  0.134 -0.425 -0.201
## convergence code: 0
## Model failed to converge with max|grad| = 0.00514629 (tol = 0.002, component 1)
 # longer closing duration in final position

2.2.5 Plateau duration di

 di.platdur1<-lmer(platdur ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = di.art)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00407375 (tol = 0.002, component 1)
 qqPlot(resid(di.platdur1))
## 9612 1135 
##  343  245
 abline(h = cut <- 30)

 di.art.plat1<-droplevels(subset(di.art, resid(di.platdur1)  < cut))
 di.platdur2<-lmer(platdur ~ position  + tenseness + (1+position|subject) + (1|word), data = di.art.plat1)
 qqPlot(resid(di.platdur2))

## 1129  835 
##  235  178
 summary(di.platdur2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: platdur ~ position + tenseness + (1 + position | subject) + (1 |  
##     word)
##    Data: di.art.plat1
## 
## REML criterion at convergence: 3029.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.68392 -0.66897 -0.01788  0.73976  2.77799 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr
##  subject  (Intercept)           34.74    5.894       
##           positionphrase-final  71.10    8.432   0.95
##  word     (Intercept)            1.87    1.368       
##  Residual                      155.79   12.482       
## Number of obs: 381, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                      Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)            34.799      2.371  8.888  14.676 1.56e-07 ***
## positionphrase-final    6.152      2.958  9.036   2.080   0.0672 .  
## tensenesslax           -2.521      1.873  1.992  -1.346   0.3111    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn  0.556       
## tensenesslx -0.397  0.004
 di.platdur.int<-lmer(platdur ~ position  * tenseness + (1+position|subject) + (1|word), data = di.art.plat1)
 anova(di.platdur2,di.platdur.int) # not better
## refitting model(s) with ML (instead of REML)
## Data: di.art.plat1
## Models:
## di.platdur2: platdur ~ position + tenseness + (1 + position | subject) + (1 | 
## di.platdur2:     word)
## di.platdur.int: platdur ~ position * tenseness + (1 + position | subject) + (1 | 
## di.platdur.int:     word)
##                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## di.platdur2     8 3054.7 3086.2 -1519.3   3038.7                         
## di.platdur.int  9 3056.2 3091.7 -1519.1   3038.2 0.4336      1     0.5102
 summary(di.platdur.int)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: platdur ~ position * tenseness + (1 + position | subject) + (1 |  
##     word)
##    Data: di.art.plat1
## 
## REML criterion at convergence: 3024.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.72167 -0.67304  0.01466  0.71742  2.80798 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr
##  subject  (Intercept)           34.763   5.896       
##           positionphrase-final  71.253   8.441   0.95
##  word     (Intercept)            1.892   1.375       
##  Residual                      156.008  12.490       
## Number of obs: 381, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                                   Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                         35.218      2.456  10.158  14.341 4.52e-08
## positionphrase-final                 5.299      3.226  12.711   1.643    0.125
## tensenesslax                        -3.351      2.254   4.114  -1.487    0.209
## positionphrase-final:tensenesslax    1.711      2.564 357.833   0.667    0.505
##                                      
## (Intercept)                       ***
## positionphrase-final                 
## tensenesslax                         
## positionphrase-final:tensenesslax    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- tnsnss
## pstnphrs-fn  0.391              
## tensenesslx -0.462  0.222       
## pstnphrs-f:  0.256 -0.397 -0.552
 # final has longer plateau phase than medial

2.2.6 Closing displacement di

 di.closdis1<-lmer(closdisp ~ position  + tenseness + (1+position|subject) + (1|word), data = di.art)
 qqPlot(resid(di.closdis1))

## 671 276 
## 152  46
 summary(di.closdis1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdisp ~ position + tenseness + (1 + position | subject) +  
##     (1 | word)
##    Data: di.art
## 
## REML criterion at convergence: 1439.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0103 -0.7152 -0.0486  0.5136  3.9933 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr 
##  subject  (Intercept)          1.2143   1.1020        
##           positionphrase-final 0.4043   0.6358   -0.29
##  word     (Intercept)          1.0049   1.0024        
##  Residual                      2.0090   1.4174        
## Number of obs: 393, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                      Estimate Std. Error      df t value Pr(>|t|)   
## (Intercept)            5.3029     0.7995  3.0638   6.632  0.00654 **
## positionphrase-final  -0.4030     0.2468  8.9582  -1.633  0.13705   
## tensenesslax           2.4667     1.0126  1.9996   2.436  0.13520   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.154       
## tensenesslx -0.633  0.000
 plot(allEffects(di.closdis1))

 di.closdis.int<-lmer(closdisp ~ position  * tenseness +
                        (1+position+tenseness|subject) + (1|word), data = di.art)
  anova(di.closdis1,di.closdis.int) # interaction for displacement better, but only slightly
## refitting model(s) with ML (instead of REML)
## Data: di.art
## Models:
## di.closdis1: closdisp ~ position + tenseness + (1 + position | subject) + 
## di.closdis1:     (1 | word)
## di.closdis.int: closdisp ~ position * tenseness + (1 + position + tenseness | 
## di.closdis.int:     subject) + (1 | word)
##                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## di.closdis1     8 1456.9 1488.7 -720.45   1440.9                             
## di.closdis.int 12 1419.5 1467.2 -697.74   1395.5 45.407      4  3.271e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# warrants closer inspection for differences in vowel height
# final is less (!) displaced than medial, other than for mono-words

2.2.7 Closing velocity di

 di.closvel1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +  (1+position+tenseness|subject) + (1|word), data = di.art)
 qqPlot(resid(di.closvel1))

## 1355  301 
##  313   61
 summary(di.closvel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + tenseness + (1 + position + tenseness |  
##     subject) + (1 | word)
##    Data: di.art
## 
## REML criterion at convergence: 2175.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2309 -0.5824 -0.0087  0.5525  2.8584 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr       
##  subject  (Intercept)          30.052   5.4820              
##           positionphrase-final  2.995   1.7306   -0.35      
##           tensenesslax          0.305   0.5523    0.03  0.18
##  word     (Intercept)           3.514   1.8744              
##  Residual                      12.808   3.5788              
## Number of obs: 393, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                      Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)           16.1868     2.2046  8.6505   7.342 5.37e-05 ***
## positionphrase-final  -2.2481     0.6558  8.9651  -3.428  0.00758 ** 
## tensenesslax           6.3497     1.9170  2.0321   3.312  0.07860 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.275       
## tensenesslx -0.429  0.014
 plot(allEffects(di.closvel1))

 di.closvel.int<-lmer(PVEL..V.cm.sec. ~ position  * tenseness + (1|subject) + (1|word), data = di.art)
 summary(di.closvel.int)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position * tenseness + (1 | subject) + (1 |  
##     word)
##    Data: di.art
## 
## REML criterion at convergence: 2175.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1417 -0.6210 -0.0223  0.5449  2.8800 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 27.673   5.261   
##  word     (Intercept)  3.487   1.867   
##  Residual             13.309   3.648   
## Number of obs: 393, groups:  subject, 10; word, 4
## 
## Fixed effects:
##                                   Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                        15.6356     2.1556   8.5889   7.254 6.11e-05
## positionphrase-final               -1.1541     0.5230 378.0212  -2.207  0.02794
## tensenesslax                        7.4360     1.9380   2.1491   3.837  0.05499
## positionphrase-final:tensenesslax  -2.1586     0.7371 378.0279  -2.929  0.00361
##                                      
## (Intercept)                       ***
## positionphrase-final              *  
## tensenesslax                      .  
## positionphrase-final:tensenesslax ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp- tnsnss
## pstnphrs-fn -0.121              
## tensenesslx -0.450  0.134       
## pstnphrs-f:  0.086 -0.710 -0.189
 anova(di.closvel1,di.closvel.int) # not better
## refitting model(s) with ML (instead of REML)
## Data: di.art
## Models:
## di.closvel.int: PVEL..V.cm.sec. ~ position * tenseness + (1 | subject) + (1 | 
## di.closvel.int:     word)
## di.closvel1: PVEL..V.cm.sec. ~ position + tenseness + (1 + position + tenseness | 
## di.closvel1:     subject) + (1 | word)
##                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## di.closvel.int  7 2196.5 2224.3 -1091.2   2182.5                         
## di.closvel1    11 2204.4 2248.1 -1091.2   2182.4 0.1016      4     0.9988
 plot(allEffects(di.closvel.int))

 di.closvel.tenseness.med<-lmer(PVEL..V.cm.sec. ~ tenseness + (1|subject) + (1|word), data = di.art %>% filter(position=="phrase-medial"))
 summary(di.closvel.tenseness.med)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ tenseness + (1 | subject) + (1 | word)
##    Data: di.art %>% filter(position == "phrase-medial")
## 
## REML criterion at convergence: 1098.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.64032 -0.60479 -0.03031  0.65028  2.84462 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 30.209   5.496   
##  word     (Intercept)  3.511   1.874   
##  Residual             12.370   3.517   
## Number of obs: 198, groups:  subject, 10; word, 4
## 
## Fixed effects:
##              Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)    15.648      2.214  8.502   7.067 7.81e-05 ***
## tensenesslax    7.420      1.939  1.998   3.826   0.0621 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tensenesslx -0.438
 di.closvel.tenseness.fin<-lmer(PVEL..V.cm.sec. ~ tenseness + (1|subject) + (1|word), data = di.art %>% filter(position=="phrase-final"))
 summary(di.closvel.tenseness.fin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ tenseness + (1 | subject) + (1 | word)
##    Data: di.art %>% filter(position == "phrase-final")
## 
## REML criterion at convergence: 1090.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9937 -0.5744 -0.0512  0.5089  2.7097 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 26.613   5.159   
##  word     (Intercept)  3.162   1.778   
##  Residual             12.981   3.603   
## Number of obs: 195, groups:  subject, 10; word, 4
## 
## Fixed effects:
##              Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)    14.493      2.092  8.331   6.927 9.98e-05 ***
## tensenesslax    5.246      1.852  2.000   2.833    0.105    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tensenesslx -0.443
 di.closvel.position.tense<-lmer(PVEL..V.cm.sec. ~ position + (1|subject) + (1|word), data = di.art %>% filter(tenseness=="tense"))
 summary(di.closvel.position.tense)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + (1 | subject) + (1 | word)
##    Data: di.art %>% filter(tenseness == "tense")
## 
## REML criterion at convergence: 1040.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.72554 -0.55207 -0.08212  0.54756  2.42999 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 27.570   5.251   
##  word     (Intercept)  6.414   2.533   
##  Residual              9.933   3.152   
## Number of obs: 195, groups:  subject, 10; word, 2
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)           15.6390     2.4628   3.2021   6.350  0.00648 **
## positionphrase-final  -1.1492     0.4523 183.0282  -2.541  0.01188 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pstnphrs-fn -0.091
 di.closvel.position.lax<-lmer(PVEL..V.cm.sec. ~ position + (1|subject) + (1|word), data = di.art %>% filter(tenseness=="lax"))
 summary(di.closvel.position.lax)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + (1 | subject) + (1 | word)
##    Data: di.art %>% filter(tenseness == "lax")
## 
## REML criterion at convergence: 1148.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6941 -0.6793  0.0318  0.5785  2.6936 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 27.8615  5.278   
##  word     (Intercept)  0.5344  0.731   
##  Residual             16.5143  4.064   
## Number of obs: 198, groups:  subject, 10; word, 2
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)           23.0729     1.7941   9.9770  12.861 1.55e-07 ***
## positionphrase-final  -3.3084     0.5781 186.0504  -5.723 4.13e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pstnphrs-fn -0.160
# final is sign. slower than medial

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

closdis.bahn<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                  data = droplevels(subset(mono.art, wordpair == "Bahn/Bann")))
 qqPlot(resid(closdis.bahn))

## 545 896 
##  68 105
plot(allEffects(closdis.bahn))

summary(closdis.bahn)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdisp ~ position + tenseness + (1 + position | subject)
##    Data: droplevels(subset(mono.art, wordpair == "Bahn/Bann"))
## 
## REML criterion at convergence: 839.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.69602 -0.57113 -0.03963  0.54048  2.17252 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr 
##  subject  (Intercept)          16.27    4.034         
##           positionphrase-final  3.97    1.993    -0.54
##  Residual                       3.02    1.738         
## Number of obs: 198, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)           16.0963     1.2934   9.1508  12.445 4.84e-07 ***
## positionphrase-final   1.1871     0.6771   9.0038   1.753    0.113    
## tensenesslax          -1.6133     0.2475 177.1990  -6.520 7.10e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.534       
## tensenesslx -0.093 -0.008
# for Bahn/Bann no effect of position, but lax vowels are less displaced
closvel.bahn<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1+position|subject),
                   data = droplevels(subset(mono.art, wordpair == "Bahn/Bann")))
qqPlot(resid(closvel.bahn))

## 11100  2013 
##    11   173
plot(allEffects(closvel.bahn))

summary(closvel.bahn)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + tenseness + (1 + position | subject)
##    Data: droplevels(subset(mono.art, wordpair == "Bahn/Bann"))
## 
## REML criterion at convergence: 1276.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.08275 -0.58251 -0.03407  0.56592  2.70483 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr 
##  subject  (Intercept)          87.14    9.335         
##           positionphrase-final 19.46    4.411    -0.35
##  Residual                      29.66    5.446         
## Number of obs: 198, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)           38.1391     3.0267   9.2941  12.601 3.75e-07 ***
## positionphrase-final  -4.7406     1.5966   8.9355  -2.969   0.0158 *  
## tensenesslax           0.4271     0.7754 177.2831   0.551   0.5824    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.361       
## tensenesslx -0.125 -0.010
# velocity is slower in final position, no effect for laxness


vh.beet<-droplevels(subset(mono.art, wordpair == "Beet/Bett"))
closdis.beet<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                   data = vh.beet)
qqPlot(resid(closdis.beet))
## 217 214 
##  27  24
abline(h = cut <- -4)

vh.beet1<-droplevels(subset(vh.beet, resid(closdis.beet) > cut))
closdis.beet1<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                    data = vh.beet1)
qqPlot(resid(closdis.beet1))

##  220 3610 
##   28    6
plot(allEffects(closdis.beet1))

summary(closdis.beet1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdisp ~ position + tenseness + (1 + position | subject)
##    Data: vh.beet1
## 
## REML criterion at convergence: 714.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.14413 -0.54919 -0.01068  0.50285  2.79077 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr 
##  subject  (Intercept)          1.9743   1.4051        
##           positionphrase-final 0.5513   0.7425   -0.70
##  Residual                      2.0013   1.4147        
## Number of obs: 193, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)            7.8967     0.4792   9.9282  16.478 1.54e-08 ***
## positionphrase-final  -0.2502     0.3111   8.8695  -0.804    0.442    
## tensenesslax           0.8208     0.2039 172.4236   4.025 8.53e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.638       
## tensenesslx -0.217  0.010
# for Beet/Bett no effect of position, but lax vowels are more displaced
closvel.beet1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1|subject),
                    data = vh.beet1)
qqPlot(resid(closvel.beet1))

## 1282 3610 
##  141    6
plot(allEffects(closvel.beet1))

summary(closvel.beet1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + tenseness + (1 | subject)
##    Data: vh.beet1
## 
## REML criterion at convergence: 1114.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2911 -0.5333 -0.0186  0.5657  3.0491 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 25.56    5.056   
##  Residual             16.39    4.048   
## Number of obs: 193, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)           22.0121     1.6792  10.2393  13.109 9.94e-08 ***
## positionphrase-final  -3.9798     0.5836 181.0077  -6.819 1.32e-10 ***
## tensenesslax           5.5692     0.5833 181.0037   9.548  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.181       
## tensenesslx -0.177  0.015
closvel.beet2<-lmer(PVEL..V.cm.sec. ~ position  * tenseness +(1|subject),
                    data = vh.beet1)
anova(closvel.beet1,closvel.beet2) # no sign. interaction
## refitting model(s) with ML (instead of REML)
## Data: vh.beet1
## Models:
## closvel.beet1: PVEL..V.cm.sec. ~ position + tenseness + (1 | subject)
## closvel.beet2: PVEL..V.cm.sec. ~ position * tenseness + (1 | subject)
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## closvel.beet1  5 1128.8 1145.1 -559.38   1118.8                         
## closvel.beet2  6 1130.0 1149.6 -559.02   1118.0 0.7139      1     0.3981
# final position is slower; lax vowels lead to faster velocity (consistent with greater displacement)


vh.ruhm<-droplevels(subset(mono.art, wordpair == "Ruhm/Rum"))
closdis.ruhm<-lmer(closdisp ~ position  + tenseness +(1|subject),
                   data = vh.ruhm)
qqPlot(resid(closdis.ruhm))

## 1671 1651 
##   74   72
plot(allEffects(closdis.ruhm))

summary(closdis.ruhm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdisp ~ position + tenseness + (1 | subject)
##    Data: vh.ruhm
## 
## REML criterion at convergence: 565.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.20549 -0.63934 -0.06795  0.54954  2.82855 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 1.1325   1.0642  
##  Residual             0.9215   0.9599  
## Number of obs: 193, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)            2.6542     0.3570  10.5163   7.435 1.68e-05 ***
## positionphrase-final  -0.1914     0.1385 181.0856  -1.382    0.169    
## tensenesslax           1.6455     0.1383 181.0378  11.902  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.189       
## tensenesslx -0.194 -0.006
# for Ruhm/Rum no effect of position, but lax vowels are more displaced
closvel.ruhm<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1|subject),
                   data = vh.ruhm)
qqPlot(resid(closvel.ruhm))
## 1671 1651 
##   74   72
abline(h = cut <- 6)

vh.ruhm.vel<-droplevels(subset(vh.ruhm, resid(closvel.ruhm) < cut))
closvel.ruhm1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1|subject),
                   data = vh.ruhm.vel)
plot(allEffects(closvel.ruhm1))

summary(closvel.ruhm1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + tenseness + (1 | subject)
##    Data: vh.ruhm.vel
## 
## REML criterion at convergence: 840.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.79612 -0.70610  0.02063  0.56983  3.13554 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 5.983    2.446   
##  Residual             4.475    2.115   
## Number of obs: 187, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)            6.4081     0.8180  10.4252   7.834 1.10e-05 ***
## positionphrase-final  -2.0843     0.3115 175.2438  -6.692 2.87e-10 ***
## tensenesslax           3.6482     0.3110 175.1901  11.732  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.188       
## tensenesslx -0.181 -0.025
closvel.ruhm1.int<-lmer(PVEL..V.cm.sec. ~ position  * tenseness +(1|subject),
                    data = vh.ruhm.vel)
anova(closvel.ruhm1,closvel.ruhm1.int) # interaction not better
## refitting model(s) with ML (instead of REML)
## Data: vh.ruhm.vel
## Models:
## closvel.ruhm1: PVEL..V.cm.sec. ~ position + tenseness + (1 | subject)
## closvel.ruhm1.int: PVEL..V.cm.sec. ~ position * tenseness + (1 | subject)
##                   Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## closvel.ruhm1      5 850.36 866.52 -420.18   840.36                           
## closvel.ruhm1.int  6 849.23 868.62 -418.62   837.23 3.1301      1    0.07686 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
vh.stiel<-droplevels(subset(mono.art, wordpair == "Stiel/still"))
closdis.stiel<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                   data = vh.stiel)
qqPlot(resid(closdis.stiel))
## 1384 1030 
##  153  115
abline(h = cut <- 5)

vh.stiel1<-droplevels(subset(vh.stiel, resid(closdis.stiel) < cut))
closdis.stiel1<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                    data = vh.stiel1)
plot(allEffects(closdis.stiel1))

summary(closdis.stiel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdisp ~ position + tenseness + (1 + position | subject)
##    Data: vh.stiel1
## 
## REML criterion at convergence: 696.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6815 -0.5611 -0.0747  0.5507  2.3478 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr 
##  subject  (Intercept)          0.9549   0.9772        
##           positionphrase-final 0.3964   0.6296   -0.35
##  Residual                      1.7865   1.3366        
## Number of obs: 195, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)            4.1797     0.3511  10.4941  11.903  2.0e-07 ***
## positionphrase-final   1.6091     0.2764   8.7510   5.822 0.000281 ***
## tensenesslax          -1.3174     0.1918 174.3271  -6.868  1.1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.415       
## tensenesslx -0.273  0.010
# for Stiel/still more displacement in final position and less displacement after lax vowels
closvel.stiel<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1+position|subject),
                    data = vh.stiel)
qqPlot(resid(closvel.stiel))
## 1024 1384 
##  109  153
abline(h = cut <- 10)

vh.stiel.vel1<-droplevels(subset(vh.stiel, resid(closvel.stiel) < cut))
closvel.stiel1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1|subject),
                     data = vh.stiel.vel1)
qqPlot(resid(closvel.stiel1))

##   612 11012 
##    71   156
plot(allEffects(closvel.stiel1))

summary(closvel.stiel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + tenseness + (1 | subject)
##    Data: vh.stiel.vel1
## 
## REML criterion at convergence: 1081.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.60653 -0.61818 -0.07522  0.54715  2.80316 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept)  9.627   3.103   
##  Residual             15.249   3.905   
## Number of obs: 191, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)            14.321      1.099  12.063  13.027 1.82e-08 ***
## positionphrase-final   -1.009      0.566 179.066  -1.783   0.0764 .  
## tensenesslax           -3.466      0.566 179.077  -6.123 5.65e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.265       
## tensenesslx -0.263  0.015

2.4 Table 6

# part 1
x<-capture.output(
  texreg::texreg(list(texreg::extract(closdis.bahn, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                      texreg::extract(closdis.beet1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                      texreg::extract(closdis.ruhm, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                      texreg::extract(closdis.stiel1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
  ),
  digits = 1,booktabs = T,custom.model.names = c("Bahn/Bann","Beet/Bett","Ruhm/Rum","Stiel/still"),
  float.pos = "!htbp", label = "tab:dispvel.mono", use.packages = F,single.row = T,caption.above = T,
  caption="LME model predictions for the articulatory measurements of closing gesture displacement (ClosDisp) and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for   phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax) in all monosyllabic target words.",
  dcolumn = T, fontsize = "footnotesize",center=F)
)
x[3]<-paste0(x[3], "\\centering")
x[15]<-paste0(x[15],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(closdis.bahn)[1],2),"/",round(r.squaredGLMM(closdis.bahn)[2],2),
              "&", round(r.squaredGLMM(closdis.beet1)[1],2),"/", round(r.squaredGLMM(closdis.beet1)[2],2),
              "&", round(r.squaredGLMM(closdis.ruhm)[1],2), "/", round(r.squaredGLMM(closdis.ruhm)[2],2),
              "&", round(r.squaredGLMM(closdis.stiel1)[1],2), "/", round(r.squaredGLMM(closdis.stiel1)[2],2),
              "\\\\ \n")
cat(x, sep = "\n")
## 
## \begin{table}[!htbp]
## \caption{LME model predictions for the articulatory measurements of closing gesture displacement (ClosDisp) and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for   phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax) in all monosyllabic target words.}\centering
## \begin{footnotesize}
## \begin{tabular}{l D{)}{)}{9)3} D{)}{)}{9)3} D{)}{)}{9)3} D{)}{)}{9)3} }
## \toprule
##  & \multicolumn{1}{c}{Bahn/Bann} & \multicolumn{1}{c}{Beet/Bett} & \multicolumn{1}{c}{Ruhm/Rum} & \multicolumn{1}{c}{Stiel/still} \\
## \midrule
## (Intercept)          & 16.1 \; (1.3)^{***} & 7.9 \; (0.5)^{***} & 2.7 \; (0.4)^{***} & 4.2 \; (0.4)^{***}  \\
## positionphrase-final & 1.2 \; (0.7)        & -0.3 \; (0.3)      & -0.2 \; (0.1)      & 1.6 \; (0.3)^{***}  \\
## tensenesslax         & -1.6 \; (0.2)^{***} & 0.8 \; (0.2)^{***} & 1.6 \; (0.1)^{***} & -1.3 \; (0.2)^{***} \\
## \midrule
## AIC                  & 853.2               & 728.1              & 575.5              & 710.8               \\
## Num. obs.            & 198                 & 193                & 193                & 195                 \\
## Num. groups: subject & 10                  & 10                 & 10                 & 10                  \\\midrule R$_m^2$/R$_c^2$ & 0.06/0.83&0.05/0.46&0.25/0.66&0.29/0.53\\ 
## 
## \bottomrule
## \multicolumn{5}{l}{\tiny{$^{***}p<0.001$, $^{**}p<0.01$, $^*p<0.05$}}
## \end{tabular}
## \end{footnotesize}
## \label{tab:dispvel.mono}
## \end{table}
# part 2
  x<-capture.output(
    texreg::texreg(list(texreg::extract(closvel.bahn, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                        texreg::extract(closvel.beet1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                        texreg::extract(closvel.ruhm1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                        texreg::extract(closvel.stiel1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
    ),
    digits = 1,booktabs = T,custom.model.names = c("Bahn/Bann","Beet/Bett","Ruhm/Rum","Stiel/still"),
    float.pos = "!htbp", label = "tab:dispvel.mono", use.packages = F,single.row = T,caption.above = T,
    caption="LME model predictions for the articulatory measurements of closing gesture displacement (ClosDisp) and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for   phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax) in all monosyllabic target words.",
    dcolumn = T, fontsize = "footnotesize",center=F)
  )
  x[3]<-paste0(x[3], "\\centering")
  x[15]<-paste0(x[15],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(closvel.bahn)[1],2),"/",round(r.squaredGLMM(closvel.bahn)[2],2),
                "&", round(r.squaredGLMM(closvel.beet1)[1],2),"/", round(r.squaredGLMM(closvel.beet1)[2],2),
                "&", round(r.squaredGLMM(closvel.ruhm)[1],2), "/", round(r.squaredGLMM(closvel.ruhm)[2],2),
                "&", round(r.squaredGLMM(closvel.stiel1)[1],2), "/", round(r.squaredGLMM(closvel.stiel1)[2],2),
                "\\\\ \n")
  cat(x, sep = "\n")
## 
## \begin{table}[!htbp]
## \caption{LME model predictions for the articulatory measurements of closing gesture displacement (ClosDisp) and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for   phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax) in all monosyllabic target words.}\centering
## \begin{footnotesize}
## \begin{tabular}{l D{)}{)}{9)3} D{)}{)}{9)3} D{)}{)}{9)3} D{)}{)}{9)3} }
## \toprule
##  & \multicolumn{1}{c}{Bahn/Bann} & \multicolumn{1}{c}{Beet/Bett} & \multicolumn{1}{c}{Ruhm/Rum} & \multicolumn{1}{c}{Stiel/still} \\
## \midrule
## (Intercept)          & 38.1 \; (3.0)^{***} & 22.0 \; (1.7)^{***} & 6.4 \; (0.8)^{***}  & 14.3 \; (1.1)^{***} \\
## positionphrase-final & -4.7 \; (1.6)^{**}  & -4.0 \; (0.6)^{***} & -2.1 \; (0.3)^{***} & -1.0 \; (0.6)       \\
## tensenesslax         & 0.4 \; (0.8)        & 5.6 \; (0.6)^{***}  & 3.6 \; (0.3)^{***}  & -3.5 \; (0.6)^{***} \\
## \midrule
## AIC                  & 1290.1              & 1124.5              & 850.1               & 1091.8              \\
## Num. obs.            & 198                 & 193                 & 187                 & 191                 \\
## Num. groups: subject & 10                  & 10                  & 10                  & 10                  \\\midrule R$_m^2$/R$_c^2$ & 0.05/0.75&0.22/0.7&0.28/0.64&0.12/0.46\\ 
## 
## \bottomrule
## \multicolumn{5}{l}{\tiny{$^{***}p<0.001$, $^{**}p<0.01$, $^*p<0.05$}}
## \end{tabular}
## \end{footnotesize}
## \label{tab:dispvel.mono}
## \end{table}
vh.di.miete<-droplevels(subset(di.art, wordpair=="Miete/Mitte"))
disp.miete<-lmer(closdisp ~ position  + tenseness + (1|subject), data = vh.di.miete)
qqPlot(resid(disp.miete))
## 671 276 
##  73  26
abline(h = cut <- 4)

vh.di.miete.sub<-subset(vh.di.miete, resid(disp.miete) < cut)
disp.miete1<-lmer(closdisp ~ position  + tenseness + (1|subject), data = vh.di.miete.sub)
qqPlot(resid(disp.miete1))

## 1350 1163 
##  146  121
summary(disp.miete1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdisp ~ position + tenseness + (1 | subject)
##    Data: vh.di.miete.sub
## 
## REML criterion at convergence: 664.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.04833 -0.55399 -0.09319  0.54075  2.71563 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 0.806    0.8978  
##  Residual             1.648    1.2837  
## Number of obs: 192, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)            4.3443     0.3257  12.7714  13.339 7.17e-09 ***
## positionphrase-final  -0.3722     0.1855 180.1272  -2.006   0.0463 *  
## tensenesslax           2.8328     0.1857 180.2153  15.251  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.275       
## tensenesslx -0.279 -0.033
# marginal effect for position; lax has greater displacement
vel.miete<-lmer(PVEL..V.cm.sec. ~ position  + tenseness + (1|subject), data = vh.di.miete)
qqPlot(resid(vel.miete))

## 1355  290 
##  154   35
plot(allEffects(vel.miete))

summary(vel.miete)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + tenseness + (1 | subject)
##    Data: vh.di.miete
## 
## REML criterion at convergence: 1093
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5672 -0.6859 -0.0171  0.6101  2.0992 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 17.54    4.188   
##  Residual             13.95    3.735   
## Number of obs: 195, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)           14.3509     1.4031  10.4544  10.228 9.00e-07 ***
## positionphrase-final  -2.2056     0.5350 183.0019  -4.122 5.68e-05 ***
## tensenesslax           7.5789     0.5353 183.0159  14.157  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.187       
## tensenesslx -0.191 -0.016
# slower velocity in final position; faster velocity after lax vowels


vh.di.huete<-droplevels(subset(di.art, wordpair=="Huete/Huette"))
disp.huete<-lmer(closdisp ~ position  + tenseness + (1|subject), data = vh.di.huete)
qqPlot(resid(disp.huete))

## 323 625 
##  40  63
summary(disp.huete)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: closdisp ~ position + tenseness + (1 | subject)
##    Data: vh.di.huete
## 
## REML criterion at convergence: 697.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.86832 -0.67487  0.05077  0.57674  2.93343 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 1.998    1.414   
##  Residual             1.690    1.300   
## Number of obs: 198, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)            6.1695     0.4750  10.5506  12.989 7.97e-08 ***
## positionphrase-final  -0.3617     0.1848 186.0091  -1.957   0.0519 .  
## tensenesslax           2.0572     0.1848 186.0108  11.130  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.196       
## tensenesslx -0.198  0.020
# no effect for position; lax has greater displacement
vel.huete<-lmer(PVEL..V.cm.sec. ~ position  + tenseness + (1|subject), data = vh.di.huete)
qqPlot(resid(vel.huete))

##  301 9212 
##   23  161
plot(allEffects(vel.huete))

summary(vel.huete)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PVEL..V.cm.sec. ~ position + tenseness + (1 | subject)
##    Data: vh.di.huete
## 
## REML criterion at convergence: 1065.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8289 -0.5045 -0.0764  0.6045  3.3512 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 40.75    6.384   
##  Residual             10.58    3.253   
## Number of obs: 198, groups:  subject, 10
## 
## Fixed effects:
##                      Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)           18.0049     2.0584   9.4780   8.747 7.66e-06 ***
## positionphrase-final  -2.2733     0.4626 186.0034  -4.914 1.95e-06 ***
## tensenesslax           5.1551     0.4626 186.0039  11.143  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pstnp-
## pstnphrs-fn -0.113       
## tensenesslx -0.115  0.020
# slower velocity in final position; faster velocity after lax vowels

2.5 Fig 6: Articulatory effect plots

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

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

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

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

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

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

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

 (ef.m.closdur.ten.plot<-
     ef.m.closdur.ten %>%
     ggplot(aes(x=tenseness, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     scale_color_brewer(palette = "Set1")+
     xlab("Closing duration (log)")+
     ylab("")+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))

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

 (ef.m.platdur.ten.plot<-
     ef.m.platdur.ten %>%
     ggplot(aes(x=tenseness, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     scale_color_brewer(palette = "Set1")+
     xlab("")+
     ylab("Plateau duration (log)")+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))

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

 (ef.m.closdis.ten.plot<-
     ef.m.closdis.ten %>%
     ggplot(aes(x=tenseness, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     scale_color_brewer(palette = "Set1")+
     xlab("")+
     ylab("Closing displacement (mm)")+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))

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

 (ef.m.closvel.ten.plot<-
     ef.m.closvel.ten %>%
     ggplot(aes(x=tenseness, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     scale_color_brewer(palette = "Set1")+
     xlab("")+
     ylab("Closing velocity (cm/s)")+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))

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

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

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

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

 ggpubr::ggarrange(ef.m.closdur.pos.plot, ef.m.platdur.pos.plot, ef.m.closdis.pos.plot, ef.m.closvel.pos.plot,
                   # ef.m.closdur.ten.plot, ef.m.platdur.ten.plot, ef.m.closdis.ten.plot, ef.m.closvel.ten.plot,
                   ef.d.closdur.pos.plot, ef.d.platdur.pos.plot, ef.d.closdis.pos.plot, ef.d.closvel.pos.plot,
                   common.legend = T,
                   labels = c("a)","b)","c)","d)","e)","f)","g)","h)"),
                   font.label = list(size=11,face="plain"),align = "h",nrow=2,ncol=4)

 ggsave(file="mono-art-effects.pdf", width = 200, height = 130, units = "mm")

2.6 Table 5

 x<-capture.output(
   texreg::texreg(list(texreg::extract(clodur2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(platdur2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(closdis2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(closvel2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
   )
   ,
   override.pvalues = list(summary(clodur2)$coefficients[1:3,5],
                           summary(platdur2)$coefficients[1:3,5],
                           summary(closdis2)$coefficients[1:3,5],
                           summary(closvel2)$coefficients[1:3,5])
   ,
   digits = 1,booktabs = T,custom.model.names = c("ClosDur (log)","PlatDur (log)","ClosDisp","ClosVel"),
   float.pos = "!htbp", label = "tab:articresults.mono", use.packages = F,single.row = T,caption.above = T,
   caption="LME model predictions for the articulatory measurements of closing gesture duration  (ClosDur),    plateau duration (PlatDur), closing gesture displacement (ClosDisp), and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for the fixed effects  phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax).",
   dcolumn = T, fontsize = "footnotesize",center=F)
 )
## Warning in override(models, override.coef, override.se, override.pvalues, :
## P-values were provided using 'override.pval', but standard errors were not
## replaced!
 x[3]<-paste0(x[3], "\\centering")
 x[16]<-paste0(x[16],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(clodur2)[1],2),"/",round(r.squaredGLMM(clodur2)[2],2),
               "&", round(r.squaredGLMM(platdur2)[1],2),"/", round(r.squaredGLMM(platdur2)[2],2),
               "&", round(r.squaredGLMM(closdis2)[1],2), "/", round(r.squaredGLMM(closdis2)[2],2),
               "&", round(r.squaredGLMM(closvel2)[1],2), "/", round(r.squaredGLMM(closvel2)[2],2),
               "\\\\ \n")
 cat(x, sep = "\n")
## 
## \begin{table}[!htbp]
## \caption{LME model predictions for the articulatory measurements of closing gesture duration  (ClosDur),    plateau duration (PlatDur), closing gesture displacement (ClosDisp), and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for the fixed effects  phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax).}\centering
## \begin{footnotesize}
## \begin{tabular}{l D{)}{)}{9)3} D{)}{)}{8)3} D{)}{)}{9)1} D{)}{)}{9)2} }
## \toprule
##  & \multicolumn{1}{c}{ClosDur (log)} & \multicolumn{1}{c}{PlatDur (log)} & \multicolumn{1}{c}{ClosDisp} & \multicolumn{1}{c}{ClosVel} \\
## \midrule
## (Intercept)          & 3.8 \; (0.1)^{***} & 3.4 \; (0.2)^{***} & 7.6 \; (2.9)^{*} & 20.5 \; (6.6)^{*}  \\
## positionphrase-final & 0.3 \; (0.0)^{***} & 0.8 \; (0.0)^{***} & 0.6 \; (0.2)^{*} & -3.2 \; (0.9)^{**} \\
## tensenesslax         & -0.1 \; (0.1)      & 0.2 \; (0.2)       & -0.1 \; (4.0)    & 1.7 \; (9.0)       \\
## \midrule
## AIC                  & 263.1              & 1052.9             & 3542.9           & 4936.5             \\
## Num. obs.            & 760                & 748                & 783              & 782                \\
## Num. groups: subject & 10                 & 10                 & 10               & 10                 \\
## Num. groups: word    & 8                  & 8                  & 8                & 8                  \\\midrule R$_m^2$/R$_c^2$ & 0.19/0.49&0.34/0.58&0/0.88&0.02/0.87\\ 
## 
## \bottomrule
## \multicolumn{5}{l}{\tiny{$^{***}p<0.001$, $^{**}p<0.01$, $^*p<0.05$}}
## \end{tabular}
## \end{footnotesize}
## \label{tab:articresults.mono}
## \end{table}

2.7 Table 7

 x<-capture.output(
   texreg::texreg(list(texreg::extract(di.clodur2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(di.platdur2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(di.closdis1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(di.closvel1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
   ),
   override.pvalues = list(summary(di.clodur2)$coefficients[1:3,5],
                             summary(di.platdur2)$coefficients[1:3,5],
                             summary(di.closdis1)$coefficients[1:3,5],
                             summary(di.closvel1)$coefficients[1:3,5]),
   digits = 1,booktabs = T,custom.model.names = c("ClosDur","PlatDur","ClosDisp","ClosVel"),
   float.pos = "!htbp", label = "tab:articresults.di", use.packages = F,single.row = T,caption.above = T,
   caption="LME model predictions for the articulatory measurements of closing gesture duration  (ClosDur),    plateau duration (PlatDur), closing gesture displacement (ClosDisp), and closing gesture velocity (ClosVel)    of the last consonant in all disyllabic word pairs (standard error in parentheses)    for the fixed effects  phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax).",
   dcolumn = T, fontsize = "footnotesize",center=F)
 )
## Warning in override(models, override.coef, override.se, override.pvalues, :
## P-values were provided using 'override.pval', but standard errors were not
## replaced!
 x[3]<-paste0(x[3], "\\centering")
 x[16]<-paste0(x[16],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(di.clodur2)[1],2),"/",round(r.squaredGLMM(di.clodur2)[2],2),
               "&", round(r.squaredGLMM(di.platdur2)[1],2),"/", round(r.squaredGLMM(di.platdur2)[2],2),
               "&", round(r.squaredGLMM(di.closdis1)[1],2), "/", round(r.squaredGLMM(di.closdis1)[2],2),
               "&", round(r.squaredGLMM(di.closvel1)[1],2), "/", round(r.squaredGLMM(di.closvel1)[2],2),
               "\\\\ \n")
 cat(x, sep = "\n")
## 
## \begin{table}[!htbp]
## \caption{LME model predictions for the articulatory measurements of closing gesture duration  (ClosDur),    plateau duration (PlatDur), closing gesture displacement (ClosDisp), and closing gesture velocity (ClosVel)    of the last consonant in all disyllabic word pairs (standard error in parentheses)    for the fixed effects  phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax).}\centering
## \begin{footnotesize}
## \begin{tabular}{l D{)}{)}{9)3} D{)}{)}{9)3} D{)}{)}{9)2} D{)}{)}{9)3} }
## \toprule
##  & \multicolumn{1}{c}{ClosDur} & \multicolumn{1}{c}{PlatDur} & \multicolumn{1}{c}{ClosDisp} & \multicolumn{1}{c}{ClosVel} \\
## \midrule
## (Intercept)          & 47.4 \; (3.4)^{***} & 34.8 \; (2.4)^{***} & 5.3 \; (0.8)^{**} & 16.2 \; (2.2)^{***} \\
## positionphrase-final & 3.4 \; (2.0)        & 6.2 \; (3.0)        & -0.4 \; (0.2)     & -2.2 \; (0.7)^{**}  \\
## tensenesslax         & 2.5 \; (4.4)        & -2.5 \; (1.9)       & 2.5 \; (1.0)      & 6.3 \; (1.9)        \\
## \midrule
## AIC                  & 2754.3              & 3045.1              & 1455.8            & 2197.5              \\
## Num. obs.            & 373                 & 381                 & 393               & 393                 \\
## Num. groups: subject & 10                  & 10                  & 10                & 10                  \\
## Num. groups: word    & 4                   & 4                   & 4                 & 4                   \\\midrule R$_m^2$/R$_c^2$ & 0.03/0.48&0.04/0.45&0.27/0.65&0.2/0.77\\ 
## 
## \bottomrule
## \multicolumn{5}{l}{\tiny{$^{***}p<0.001$, $^{**}p<0.01$, $^*p<0.05$}}
## \end{tabular}
## \end{footnotesize}
## \label{tab:articresults.di}
## \end{table}
 # plot articulatory data in boxplots from data subsets used for the models above
 # mono.art.sub1 # clodur mono
 # mono.art.plat1 # platdur mono
 # mono.art # clodis mono
 # mono.art.vel1 # clovel mono
 # di.art.sub1 # clodur di
 # di.art.plat1 # platdur di
 # di.art # clodis di
 # di.art # clovel di

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

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


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

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

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

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


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

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

 ggpubr::ggarrange(m.clos,m.plat,m.disp,m.vel,
                   common.legend = T,
                   labels = c("a)","b)","c)","d)"),
                   font.label = list(size=11,face="plain"),align = "h",nrow=2,ncol=2)

 ggsave(file="mono-art-raw-tenseness-on-x.pdf", width = 250, height = 250, units = "mm")


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

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

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


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

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

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

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

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


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

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



 ggpubr::ggarrange(d.clos,d.plat,d.disp,d.vel,
                   common.legend = T,
                   labels = c("a)","b)","c)","d)"),
                   font.label = list(size=11,face="plain"),align = "h",nrow=2,ncol=2)

 ggsave(file="di-art-raw-tenseness-on-x.pdf", width = 200, height = 200, units = "mm")

3 Figure of proportional lengthening

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