The human voice aligns with whole-body kinetics

Supporting information

Show code setting up the packages

#R-packages
library(papaja)   #for Rmarkdown template for APA(ish) manuscript
library(distill) #for html output
library(ggplot2)    #for plotting
library(gridExtra)  #for plotting in panels
library(knitr)  #for document generation
library(magick)  #for plots
library(cowplot) #for plots
library(plyr) #for revalue
library(tinytex) #for outputting pdfs
library(readr)  #for data reading
library(bookdown) #for crossreferencing
library(rstudioapi) #for setting file paths
library(kableExtra) #for producing nicer tables
library(dplyr) #for data wrangling
library(ggbeeswarm) #for beeswarm plots
library(modelsummary) #for handling model summary tables
library(lme4) #for statistical modeling
library(lmerTest) #for statistical modeling
library(emmeans) #for post-hoc analysis
library(broom.mixed) #for turning fitted models into tidy data frames

r_refs("references.bib")

This document contains further information on the supporting methods and results of “The human voice aligns with whole-body kinetics”. For clarity, it is indicated on the right, what sections in the paper the respective part refers to. Code chunks are hidden by default, but can be made visible by clicking “Show code”.

Supporting methods and results

This is a fully computationally reproducible manuscript written in R Markdown. The full dataset is available at the Donders Repository (redacted for review).

This refers to the section “Supporting information” of the paper.

Show code loading and preparing the data

#Local data
localfolder <- "C:/Users/U941163/Desktop/Raphael/projects/VENI-local/" #THIS NEEDS TO BE CHANGED

#lets set all the data folders
rawd      <- paste0(localfolder, 'Trials/')    #raw trial level data
procd     <- paste0(localfolder, 'Processed/triallevel/') #processed data folder
procdtot  <- paste0(localfolder, 'Processed/complete_datasets/') #processed data folder
daset     <- paste0(localfolder, 'Dataset/')   #dataset folder
meta      <- paste0(localfolder, 'Meta/')   #Meta and trial data

#load in the trialinfo data (containing info like condition order, trial number etc)
trialf    <- list.files(meta, pattern = 'triallist*')
triallist <- data.frame()
for(i in trialf)
{
  tr        <- read.csv(paste0(meta, i), sep = ';')
  tr$ppn    <- parse_number(i)
  triallist <- rbind.data.frame(triallist,tr)
}

#make corrections to trial list due to experimenter errors
  # p10 trial 21 run with a weight
triallist$weight_condition[triallist$ppn==10 & triallist$trial=="21"] <- 'weight' #triallist$trialindex[triallist$ppn==10 & triallist$trial=="21"]
  # p14 double check trial "51", "52", "53" were all weight condition
triallist$weight_condition[triallist$ppn==14 & triallist$trial%in% c("51", "52", "53")] <- 'weight' #triallist$trialindex[triallist$ppn==14 & triallist$trial%in% c("51", "52", "53")]


#load in the metainfo data (containing info like gender, handedness, etc)
mtf     <- list.files(meta, pattern = 'bodymeta*')
mtlist  <- data.frame()
for(i in mtf)
{
  tr      <- read.csv(paste0(meta, i), sep = ';')
  mtlist  <- rbind.data.frame(mtlist,tr)
}
mtlist$ppn <- parse_number(as.character(mtlist$ppn))

Supporting methods

This study has been approved by the Ethics Committee Social Sciences (ECSS) of the Radboud University (reference nr.: 22N.002642).

Show code adding participant and condition information

# add some info about the dataset
participants  <- c(unique(mtlist$ppn))
numpart       <- length(participants)

# number & percentage female & male participants
part_fem <- round(sum(mtlist$sex=='f'), 1)
perc_fem  <- round((sum(mtlist$sex=='f')/numpart)*100, 1)
part_male <- round(sum(mtlist$sex=='m'), 1)
perc_male <- round((sum(mtlist$sex=='m')/numpart)*100, 1)

# body info
m_age           <- round(mean(mtlist$age), 1)
sd_age          <- round(sd(mtlist$age), 1)
m_weight        <- round(mean(mtlist$weight), 1)
sd_weight       <- round(sd(mtlist$weight), 1)
m_height        <- round(mean(mtlist$height), 1)
sd_height       <- round(sd(mtlist$height), 1)
bmi             <- mtlist$weight/((mtlist$height/100)^2)
m_bmi           <- round(mean(bmi), 1)
sd_bmi          <- round(sd(bmi), 1)

# other body measurements
  # Function to calculate the mean of a character vector
  mean_of_char_vector <- function(vec) {
    mean(as.numeric(vec))}

skinfold  <- strsplit(mtlist$upperarmfold, '_')
# Applying the function to each element of the list and storing the results in a new vector
skinfs    <- sapply(skinfold, mean_of_char_vector)
m_skin    <- round(mean(skinfs), 1)
sd_skin   <- round(sd(skinfs), 1)
# we also measure upper/under arm length (but will leave it at this)

# read data
tr_wd <- read.csv(paste0(procdtot, 'fulldata.csv'))
tr_wd <- subset(tr_wd, trialindex > 8 & vocal_condition=='vocalize') #only real trials and only vocalization
numtrials = nrow(tr_wd)
perc_w = round((sum(tr_wd$weight_condition=='weight')/numtrials)*100, 1) #percentage weight
perc_pas = round((sum(tr_wd$movement_condition=='no movement')/numtrials)*100, 1) #percentage passive
perc_int = round((sum(tr_wd$movement_condition=='internal rotation')/numtrials)*100, 1) #percentage internal rotation
perc_ext = round((sum(tr_wd$movement_condition=='external rotation')/numtrials)*100, 1) #percentage external rotation
perc_flex = round((sum(tr_wd$movement_condition=='flexion')/numtrials)*100, 1) #percentage flexion
perc_extens = round((sum(tr_wd$movement_condition=='extension')/numtrials)*100, 1) #percentage extension

Experimental design

This study concerned a two-level wrist-weight manipulation (no weight vs. weight), a two-level within-subject vocalization condition (expire vs. vocalize), and a five-level within-subject movement condition (‘no movement’, ‘extension’, ‘flexion’, ‘external rotation’, ‘internal rotation’). With 4 trial repetitions over the experiment, we yield 80 (2 weights x 2 vocalizations x 5 movements x 4 repetitions) trials per participant. Trials were blocked by weight condition and vocalization condition (so that weights and task did not switch from trial to trial). Within blocks all movement conditions were randomized.

Participants

For the current pre-registered confirmatory experiment, as planned and supported by a power analysis (see preregistration), we collected N = (17) participants: 7 female, 10 male, M (SD) age = 28.50 (6.50), M (SD) body weight = 72.10 kg (10.20), M (SD) body height = 175.10 cm (8.50), M (SD) BMI = 23.40 (2.20), M (SD) triceps skinfold = 19.10 mm (4.30).

Exclusions and deviations from pre-registration

We also performed the experiment with one other participant, but due to an issue with LSL streaming, this dataset could not be synchronized and was lost. Furthermore, due to running over time one participant had to terminate the study earlier about half way through. Note further, that in our pre-registration we wanted to admit participants with a BMI lower than 25, but since participants were difficult to recruit we accepted three participants with slightly higher BMIs too (max BMI of the current dataset: 27.10). Participants were all able-bodied and did not have any constraints in performing the task. Finally, as stated in the pre-registration we will only report results on the vocalization trials (and not the expiration-only trials). This means that for this report in total we have N = 17 participants, with 636 analyzable trials, with balanced conditions: weight = 50.20%, no movement = 19.70%, internal rotation = 20.00 %, external rotation = 20.30%, flexion = 20.10 %, extension = 20.00 %.

Measurements and equipment

Body measurements

To enable future analyses of possible factors modulating individual-specific body properties, we collect some basic information about body properties. Namely, weight, under arm length, upper arm length, triceps skinfold, and upper arm circumference.

Experiment protocol

The experiment was coded in Python using functions from PsychoPy. The experiment was controlled via a Brainvision Button Box (Brain Products GmbH, Munich, Germany), which was also streaming its output to the data collection PC unit.

Wrist weights

To manipulate the mass set in motion, we apply a wrist weight. We use a TurnTuri sports wrist weight of 1 kg.

Video and kinematics

The participants are recorded via a videocamera (Logitech StreamCam), sampling at 60 frames per second. We used Mediapipe (Lugaresi et al. 2019) to track the skeleton and facial movements, which is implemented in Masked-piper which we also use for masking the videos (Owoyele et al. 2022). The motion-tracked skeleton, specifically the wrist of the dominant hand, is used to estimate movement initiation, peak speed, and the end of the movement. The motion tracking is, however, only used for determining movement windows, and is not of central concern.

Muscle activity (Surface ElectroMyography: sEMG)

We measured sEMG using a wired BrainAmp ExG system (Brain Products GmbH, Munich, Germany). Disposable surface electrodes (Kendall 24mm Arbo H124SG) were used, and for each of the four muscle targets we had 3 (active, reference, ground) electrodes (12 electrodes total). The sEMG system sampled at 2,500 Hz (for post-processing filters see below).
For an overview of the electrode attachments, see Fig. S1. We prepare the skin surface for EMG application with a scrub gel (NuPrep) followed by cotton ball swipe with alcohol (Podior 70 %). Active and reference electrodes were attached with a 15mm distance center to center.

Figure S1: Overview sEMG target muscles. Active (a) and reference (r), and ground (g) sEMG electrodes, for each muscle target.

This refers to Figure S1 in the paper.

We attached electrodes for focal muscles which directly participate in the internal (pectoralis major) and external rotation (infraspinatus) of the humerus. Electrodes were applied for focal muscles ipsilaterally (relative to the dominant hand). We attached electrodes to the muscle belly of the clavicular head of the pectoralis major, with a ground electrode on the clavicle on the opposite side.
We also attached electrodes for postural muscles which will likely anticipate and react to postural perturbations due to upper limb movements. Since these muscles should act in the opposite direction of the postural perturbation of the dominant hand, we applied electrodes contralaterally to the dominant hand. We attach electrodes to the rectus abdominis, with a ground electrode on the iliac crest on the opposite side. We also attached electrodes to the erector spinae muscle group (specifically, the iliocostalis lumborum).

Ground reaction measurements

We used an inhouse-built 1m² balance board with vertical pressure sensors. The sensors were derived and remodified from the Wii-Balance board sensors. The sampling rate was 400 Hz. The system was time-locked within millisecond accuracy, and has a spatial accuracy of several sub-millimeters. A national instruments card, USB-62221 performed the A/D conversion and was connected via USB to the PC.

Acoustics

To ensure proper acoustic intensity measurements we used a headset microphone; MicroMic C520 (AKG, Inc.) headset condenser cardioid microphone sampling at 16 kHz. The gain levels of the condenser power source were set by the hardware (and could not be changed).

Here are five audio examples from a male participant producing the sustained vocalization in different conditions without added weight.

Show code setting up the audio examples

library(htmltools)

audioexamples <- paste0(basefolder, "/AudioExamples/")

audioexample_nomovement       <- paste0(audioexamples, "sample1_nomovement.wav")
audioexample_internalrotation <- paste0(audioexamples, "sample2_internalrotation.wav")
audioexample_flexion          <- paste0(audioexamples, "sample3_flexion.wav")
audioexample_externalrotation <- paste0(audioexamples, "sample4_externalrotation.wav")
audioexample_extension        <- paste0(audioexamples, "sample5_extension.wav")

html_tag_audio <- function(file, type = c("wav")) {
  type <- match.arg(type)
  htmltools::tags$audio(
    controls = NA,
    htmltools::tags$source(
      src = file,
      type = glue::glue("audio/{type}", type = type)
    )
  )
}

The first example is from the no movement condition:

The second example is from the internal rotation condition:

The third example is from the external rotation condition:

The fourth example is from the extension condition:

The fifth example is from the flexion condition:

Recording setup and synchronization

We use LabStreamLayer that provides a uniform interface for streaming different signals along a network, where a common time stamp for each signal ensures sub-millisecond synchronization. We used a Linux system to record and stream the microphone recordings. Additionally a second PC collected video, and streamed ground reaction forces, and EMG. A data collection PC collected the audio, ground reaction force, and EMG streams and stored the output in XDF format for efficient storing of multiple time-varying signals. The video data was synchronized by streaming the frame number to the LSL recorder, allowing us to match up individual frames with the other signals (even when a frame is dropped).

Procedure

Participants are admitted to the study based on exclusion criteria and sign an informed consent. We ask participants to take off their shoes and we proceed with the body measurements, while instructing the participant about the nature of the study. After body measurements, we apply the surface EMG. We prepared the muscle site with rubbing gel and alcohol, and active/reference electrodes were placed with a distance of 15 mm from each other’s center. See Fig. S1 for the sEMG electrode locations. The procedures up to the start of the experiment take about 20 minutes or less in total.

Upon the start of the experiment participants take a standing position on the force platform. The experiment commences with calibration and practice trials. First, 10 seconds of silent breathing with body movements are recorded. Then participants are asked to take a maximum inspiration followed by a maximum expiration to measure signal conditions under respiratory boundary conditions. Then, for the practice trials, each movement was practiced with expiring and vocalization while performing the movement conditions, and the participant is introduced to wearing the wrist weight of 1 kg. After practice trials, participants performed 80 blocked trials.

For each (practice) trial participants are closely guided by the information on the monitor. Firstly, participants are shown the movement to be performed for the trial, and have to prompt the experimenter that they are ready to continue. Then participants are instructed to adopt the start position of the movement, which is a 90 degree elbow flexion, with either an externally rotated humerus (start position for internal rotation), or a non-rotated humerus with the wrist in front of the body (rest position for the other movement conditions). For the no movement condition participants are asked to rest their arms alongside their body. Upon trial start, participants are asked to inhale deeply with a timer counting down from 4 seconds. Then, participants are asked to ‘vocalize’ or ‘expire’, with a screen appearing after 3 seconds to perform the movement with visual guidance to where the movement end position is so that participants are reminded of the movement. After an additional 4 seconds the trial ends, which allows more than enough time to perform the movement and stabilize vocalization after the perturbation. In the no movement condition, a prompt is given to maintain one’s posture.

Preprocessing of the data streams

EMG

To reduce heart rate artifacts we apply a common method (Drake and Callaghan 2006) of high-pass filtering the signal at 30 Hz using a zero-phase 4th order Butterworth filter. We then full-wave rectified the EMG signal and applied a zero-phase low-pass 4th order Butterworth filter at 20 Hz. When filtering any signal we pad the signals to avoid edge effects. We normalized the EMG signals within participants before submitting to analyses.

Show code preprocessing EMG data

# load in the data and save it
library(signal)
library(ggrepel)

#####################preprocessing functions
butter.it <- function(x, samplingrate, order, cutoff, type = "low")
{
  nyquist <- samplingrate/2
  xpad <- c(rep(0, 1000), x, rep(0,1000)) #add some padding
  bf <- butter(order,cutoff/nyquist, type=type) 
  xpad <- as.numeric(signal::filtfilt(bf, xpad))
  x <- xpad[1000:(1000+length(x)-1)] #remove the padding
}

#rectify and high and low pass EMG signals
reclow.it <- function(emgsignal, samplingrate, cutoffhigh, cutofflow)
  {
    #high pass filter
    output <- butter.it(emgsignal, samplingrate=samplingrate, order=4, cutoff = cutoffhigh,
                        type = "high")
    #rectify and low pass
    output <<- butter.it(abs(output), samplingrate=samplingrate,  order=4, cutoff= cutofflow,
                        type = "low")
}
###########################################

####################################LOAD IN EMG DATA
EMGdata <- list.files(rawd, pattern = '*Dev_1.csv')
####################################SET common colors for muscles
col_pectoralis    <- '#e7298a'
col_infraspinatus <- '#7570b3' 
col_rectus        <- '#d95f02'
col_erector       <- '#1b9e77'
colors_mus <- c("pectoralis major" = col_pectoralis, "infraspinatus" = col_infraspinatus, "rectus abdominis" = col_rectus, "erector spinae" = col_erector)

#show an example
emgd <- read.csv(paste0(rawd, EMGdata[822]))
colnames(emgd) <- c('time_s', 'infraspinatus', 'rectus_abdominis', 'pectoralis_major', 'erector_spinae', 'empty1', 'empty2', 'empty3', 'empty4', 'respiration_belt', 'button_box')
samplingrate <- 1/mean(diff(emgd$time_s))
nyquistemg <- samplingrate/2

#example without strong high pass filter
emgd$time_s <- emgd$time_s-min(emgd$time_s)
emgd$pectoralis_major_sm <- reclow.it(emgd$pectoralis_major, samplingrate= samplingrate, cutoffhigh = 2, cutofflow = 20)
emgd$infraspinatus_sm <- reclow.it(emgd$infraspinatus, samplingrate= samplingrate, cutoffhigh = 2, cutofflow = 20)
emgd$erector_spinae_sm <- reclow.it(emgd$erector_spinae, samplingrate= samplingrate, cutoffhigh = 2, cutofflow = 20)
emgd$rectus_abdominis_sm <- reclow.it(emgd$rectus_abdominis, samplingrate= samplingrate, cutoffhigh = 2, cutofflow = 20)

#example trial
unfilteredEMG <- ggplot(emgd[seq(1, nrow(emgd),by=5),], aes(x=time_s))+geom_path(aes(y=infraspinatus_sm, color = 'infraspinatus'))+
  geom_path(aes(y=pectoralis_major_sm, color = 'pectoralis major'))+
  geom_path(aes(y=rectus_abdominis_sm, color = 'rectus abdominis'))+
  geom_path(aes(y=erector_spinae_sm, color = 'erector spinae'))+ 
  geom_text(aes(x = 0.5, y = 600), color = "black", label = "Heart rate peak 1", angle=45, size = 2) +
  geom_text(aes(x = 1.4, y = 600), color = "black", label = "Heart rate peak 2", angle=45, size = 2) +
  geom_text(aes(x = 4.4, y = 600), color = "black", label = "Heart rate peak 5", angle=45, size = 2) +
  scale_color_manual(values = colors_mus)
unfilteredEMG <- unfilteredEMG +
  labs(x = "time in seconds", y = "EMG rectified") +
  xlab('time in seconds') +
  ylab('EMG rectified') +
  ylim(0, 800) +
  theme_cowplot(12) +
  theme(legend.position="none")


#example without strong high pass filter
emgd$pectoralis_major_sm_h <- reclow.it(emgd$pectoralis_major, samplingrate= samplingrate, cutoffhigh = 30, cutofflow = 20)  
emgd$infraspinatus_sm_h  <- reclow.it(emgd$infraspinatus, samplingrate= samplingrate, cutoffhigh = 30, cutofflow = 20) 
emgd$erector_spinae_sm_h  <- reclow.it(emgd$erector_spinae, samplingrate= samplingrate, cutoffhigh = 30, cutofflow = 20) 
emgd$rectus_abdominis_sm_h  <- reclow.it(emgd$rectus_abdominis, samplingrate= samplingrate, cutoffhigh = 30, cutofflow = 20) 

#example trial
filteredEMG <- ggplot(emgd[seq(1, nrow(emgd),by=5),], aes(x=time_s)) +
  geom_path(aes(y=infraspinatus_sm_h, color = 'infraspinatus')) +
  geom_path(aes(y=pectoralis_major_sm_h, color = 'pectoralis major')) +
  geom_path(aes(y=rectus_abdominis_sm_h, color = 'rectus abdominis')) +
  geom_path(aes(y=erector_spinae_sm_h, color = 'erector spinae')) +
  scale_color_manual(values = colors_mus)
filteredEMG <- filteredEMG +
  labs(x = "time in seconds", y = "EMG rectified", color = "Legend") +
  xlab('time in seconds') +
  ylab('EMG rectified') +
  theme_cowplot(12)+
  theme(
    legend.position = c(1, 1),
    legend.justification = c("right", "top"),
    legend.box.just = "right",
    legend.margin = margin(4, 4, 4, 4)
    )

Figure S2: Example of smoothing settings for EMG signals. The upper panel shows the raw rectified high-pass filtered EMG, and the lower panel shows the low-pass filtered data to reduce artifacts of heart rate. This example shows an internal rotation trial, where we successfully retrieve the peak in pectoralis major that internally rotates the arm.

This refers to Figure S2 in the paper.

Ground reaction forces

We upsampled the balanceboard from 400 Hz to 2,500 Hz. We then applied a zero-phase low-pass 20 Hz 2nd order Butterworth filter to the padded signals. As a key measure for postural perturbation we computed the change in 2D magnitude (L2 norm of the center of pressure x and y) in center of pressure (hereafter COPc).

Acoustics

For acoustics we extract the smoothed amplitude envelope (hereafter envelope). For the envelope we apply a Hilbert transform to the waveform signal, then take the complex modulus to create a 1D time series, which is then resampled at 2,500 Hz, and smoothed with a Hann filter based on a Hanning Window of 12 Hz. We normalize the amplitude envelope signals within participants before submitting to analyses.

Show code extracting acoustics

library(rPraat) #for reading in sounds
library(dplR)   #for applying Hanning
library(seewave) #for signal processing general
library(wrassp) #for acoustic processing

##################### MAIN FUNCTION TO EXTRACT SMOOTHED ENVELOPE ###############################
amplitude_envelope.extract <- function(locationsound, smoothingHz, resampledHz)
{
  #read the sound file into R
  snd <- rPraat::snd.read(locationsound)
  #apply the hilbert on the signal
  hilb <- seewave::hilbert(snd$sig, f = snd$fs, fftw =FALSE)
  #apply complex modulus
  env <- as.vector(abs(hilb))
  #smooth with a hanning window
  env_smoothed <- dplR::hanning(x= env, n = snd$fs/smoothingHz)
  #set undeterminable at beginning and end NA's to 0
  env_smoothed[is.na(env_smoothed)] <- 0
  #resample settings at desired sampling rate
  f <- approxfun(1:(snd$duration*snd$fs),env_smoothed)
  #resample apply
  downsampled <- f(seq(from=0,to=snd$duration*snd$fs,by=snd$fs/resampledHz))
  #let function return the downsampled smoothed amplitude envelope
  return(downsampled[!is.na(downsampled)])
}

Data aggregation

All signals were sampled at, or upsampled to, 2,500 Hz. Then we aggregated the data by aligning time series in a combined by-trial format to increase ease of combined analyses. We linearly interpolated signals when sample times did not align perfectly.

Show code aggregating the data

library(pracma)
library(zoo)

overwrite = FALSE

#lets loop through the triallist and extract the files
triallist$uniquetrial <- paste0(triallist$trialindex, '_', triallist$ppn) 
  #remove duplicates so that we ignore repeated trials
howmany_repeated <- sum(duplicated(triallist$uniquetrial))
triallist <- triallist[!duplicated(triallist$uniquetrial),]


for(ID in triallist$uniquetrial)
{
  #print(paste0('working on ', ID))
  #debugging
  #ID <- triallist$uniquetrial[27] #I commented that out bc it was so many lines
  
  triallist_s <- triallist[triallist$uniquetrial==ID,]
  ##load in row triallist
  pp <- triallist_s$ppn
  pps <- triallist_s$ppn
  if((pp == '41') | (pp == '42')){pps <- '4'}
  gender <- mtlist$sex[mtlist$ppn==pps]
  hand   <- mtlist$handedness[mtlist$ppn==pps]
  triali <- triallist_s$trialindex
  
  if(!file.exists(paste0(procd, 'pp_', pp, '_trial_', triali, '_aligned.csv')) | overwrite == TRUE)
  {
  ######load in EMG and resp belt
    EMGr <- read.csv(paste0(rawd, 'p', pp, '_trial_',triali, '_BrainAmpSeries-Dev_1.csv'))
    channels <- c('time_s', 'infraspinatus', 'rectus_abdominis', 'pectoralis_major', 'erector_spinae', 'empty1', 'empty2', 'empty3', 'empty4', 'respiration_belt', 'button_box')
    colnames(EMGr) <- channels
    #only keep the relevant vectors
    EMGr <- cbind.data.frame(EMGr$time_s,EMGr$respiration_belt, EMGr$infraspinatus,
                             EMGr$rectus_abdominis,EMGr$pectoralis_major, EMGr$erector_spinae)
    colnames(EMGr) <- channels[c(1, 10, 2, 3, 4, 5)] #rename
    EMGsamp <- 1/mean(diff(EMGr$time_s-min(EMGr$time_s)))
    #smooth
      #note the respiration belt is not installed now, but if we resolve technical issues we might add it
    EMGr$respiration_belt <- butter.it(EMGr$respiration_belt, samplingrate = EMGsamp, order = 2, cutoff = 30)
      #smooth EMG with described settings
    EMGr[,3:6] <- apply(EMGr[,3:6], 2, function(x) reclow.it(scale(x, scale = FALSE, center =TRUE), samplingrate= EMGsamp, 
                                                             cutoffhigh = 30, cutofflow = 20)) #CHANGECONF we center before we smooth now
  #########load in acoustics
  locsound <- paste0(rawd, 'p', pp, '_trial_',triali, '_Micdenoised.wav')
  if(!file.exists(locsound)){print(paste0('this mic file does not exist: ', locsound))}
  if(file.exists(locsound)){
    
  
  env <- amplitude_envelope.extract(locationsound=locsound, smoothingHz =  12, resampledHz = 2500)
  acoustics <- cbind.data.frame(env)
  acoustics$time_s <- seq(0, (nrow(acoustics)-1)*1/2500, 1/2500)
  duration_time <- (max(acoustics$time_s)-min(acoustics$time_s)) #note down duration time of this trial
  #########load in motion
  mt <- read.csv(paste0(rawd, 'p', pp, '_trial_',triali, '_video_raw_body.csv'))
  if(hand == 'r'){mov <- cbind.data.frame(mt$X_RIGHT_WRIST, mt$Y_RIGHT_WRIST, mt$Z_RIGHT_WRIST)}  #CHANGECONF change we also take depth
  if(hand == 'l'){mov <- cbind.data.frame(mt$X_LEFT_WRIST, mt$Y_LEFT_WRIST, mt$Z_RIGHT_WRIST)}    #CHANGECONF we also take depth
  colnames(mov) <- c('x_wrist', 'y_wrist', 'z_wrist')
  mov$y_wrist <- mov$y_wrist*-1 #flip axes so that higher up means higher values (check)
  #filter movements
  mov[,1:3] <- apply(mov[,1:3], 2, FUN=function(x) butter.it(x, order=4, samplingrate = 60, cutoff = 10, type='low'))
  
  #########load in balanceboard
  bb <- read.csv(paste0(rawd, 'p', pp, '_trial_',triali, '_BalanceBoard_stream.csv'))
  colnames(bb) <- c('time_s', 'left_back', 'right_forward', 'right_back', 'left_forward')
  bbsamp <- 1/mean(diff(bb$time_s - min(bb$time_s)))
  bb[,2:5] <- apply(bb[,2:5], 2, function(x) sgolayfilt(x, 5, 51))
  COPX  <- (bb$right_forward+bb$right_back)-(bb$left_forward+bb$left_back)
  COPY  <- (bb$right_forward+bb$left_forward)-(bb$left_back+bb$left_back)
  bb$COPXc <- sgolayfilt(c(0, diff(COPX)), 5, 51)
  bb$COPYc <-  sgolayfilt(c(0, diff(COPY)), 5, 51)
  bb$COPc <- sqrt(bb$COPXc^2 + bb$COPYc^2)
  #########combine everything
    #combined acoustics and EMG(+resp)
  EMGr$time_ss <- EMGr$time_s-min(EMGr$time_s)
  ac_emg <- merge(acoustics, EMGr, by.x = 'time_s', by.y = 'time_ss', all=TRUE)
  ac_emg[,4:9] <- apply(ac_emg[,4:9], 2, function(y) na.approx(y, x=ac_emg$time_s, na.rm = FALSE))
  ac_emg <- ac_emg[!is.na(ac_emg$env),]
  ac_emg$time_s <- ac_emg$time_s + min(ac_emg[,4],na.rm=TRUE) #get the original time
  ac_emg <- ac_emg[!is.na(ac_emg$time_s),-4] #remove the centered time variable, and remove time NA at the end
    #now add balance board
  ac_emg_bb <- merge(ac_emg, bb[,c(1, 6:8)], by.x='time_s', by.y='time_s', all= TRUE)
  ac_emg_bb[,9:11] <- apply(ac_emg_bb[,9:11], 2, function(y) na.approx(y, x=ac_emg_bb$time_s, na.rm = FALSE))
  ac_emg_bb <- ac_emg_bb[!is.na(ac_emg_bb$env), ]
  ac_emg_bb <- ac_emg_bb[!is.na(ac_emg_bb$time_s),-4]
    #now add the final motion (which is thereby upsampled to 2500 hertz)
  start_time <- min(ac_emg_bb$time_s,na.rm=TRUE)
  mov$time_s <-   start_time+seq(0, duration_time-(duration_time /nrow(mov)), by= duration_time /nrow(mov))
  ac_emg_bb_m <- merge(ac_emg_bb, mov, by.x='time_s', by.y='time_s', all = TRUE)
  ac_emg_bb_m[,11:13] <- apply(ac_emg_bb_m[,11:13], 2, function(y) na.approx(y, x=ac_emg_bb_m$time_s, na.rm = FALSE))
  ac_emg_bb_m <- ac_emg_bb_m[!is.na(ac_emg_bb_m$env), ]
  #now write the files away to a processed folder
  write.csv(ac_emg_bb_m, paste0(procd, 'pp_', pp, '_trial_', triali, '_aligned.csv'))
    }
  }
}

Data sharing & Privacy

Video data is deidentified using the masked-piper tool to mask faces and body while maintaining kinematic information (Owoyele et al. 2022).

Overview data analyses

Note that there is a general decline in the amplitude of the vocalization during a trial (see Fig. S3, top panel). This is to be expected, as the subglottal pressure falls when the lungs deflate. To quantify deviations from stable vocalizations, we therefore detrend the amplitude envelope time series, so as to assess positive or negative peaks relative to this trend line. For the envelope, muscle activity, and the change in center of pressure we will measure the global maxima happening within the analyses window (i.e., within a trial we take a local maximum occurring between movement onset and offset). We will analyze positive and negative peaks within the movement window separately.

Show code scaling data within participant

#We have to do some within participant scaling, such that
  #all EMG's are rescaled within participants
fs <-  list.files(procd)
fs <- fs[!fs%in%list.files(procd, pattern = 'zscal_*')]

#set to true if you want to overwrite files (to check [modified] code)
overwrite = FALSE


if((length(list.files(procd, pattern = 'zscal_*')) == 0) | (overwrite==TRUE))
{
  mr <- data.frame()
  for(f in fs)
  {
    temp <- read.csv(paste0(procd, f))
    temp$pp <- parse_number(f)
    temp$tr <- f
    mr <- rbind.data.frame(mr, temp)
  }
  
  #set ppn to 4
  temp$pp[temp$pp==41 | temp$pp==42] <- 4
  
  #center the muscle measurements for each trial (due to drift)
  mr$infraspinatus <- ave(mr$infraspinatus, mr$tr, FUN = function(x){x<-x-mean(x, na.rm=TRUE)})       
  mr$rectus_abdominis <- ave(mr$rectus_abdominis, mr$tr, FUN = function(x){x<-x-mean(x, na.rm=TRUE)}) 
  mr$pectoralis_major <- ave(mr$pectoralis_major, mr$tr, FUN = function(x){x<-x-mean(x, na.rm=TRUE)}) 
  mr$erector_spinae <- ave(mr$erector_spinae, mr$tr, FUN = function(x){x<-x-mean(x, na.rm=TRUE)})     

  
  #normalize signals by participant
  mr$infraspinatus <- ave(mr$infraspinatus, mr$pp, FUN = function(x)scale(x, center = FALSE))       
  mr$rectus_abdominis <- ave(mr$rectus_abdominis, mr$pp, FUN = function(x)scale(x, center = FALSE)) 
  mr$pectoralis_major <- ave(mr$pectoralis_major, mr$pp, FUN = function(x)scale(x, center = FALSE)) 
  mr$erector_spinae <- ave(mr$erector_spinae, mr$pp, FUN = function(x)scale(x, center = FALSE))     
  mr$env_z <- ave(mr$env, mr$pp, FUN = function(x)scale(x, center = FALSE))
  
  #save the scaled variables in new objects
  for(f in unique(mr$tr))
  {
    sb <- mr[mr$tr==f, ] #subset
    #center motion per trial
    sb$x_wrist <- sb$x_wrist-mean(sb$x_wrist,na.rm=TRUE)
    sb$y_wrist <- sb$y_wrist-mean(sb$y_wrist,na.rm=TRUE)
    sp <- c(0, sqrt(diff(sb$x_wrist)^2+diff(sb$y_wrist)^2))
    sp[is.na(sp)] <- 0 #there are some trailing NA's which we fill with 0 
    sb$speed <- butter.it(sp, samplingrate = 2500, order = 4, cutoff = 5)
    sb$time_s   <- sb$time_s #this is the time
    sb$time_s_c <- sb$time_s-min(sb$time_s)# accumulating at trial start (resyncing) resynced with the psychopy software
    write.csv(sb, paste0(procd, 'zscal_', f))
  }
}

Show code producing the plot

ptest=8
ex2 <- triallist[triallist$vocal_condition=='vocalize' & triallist$movement_condition=='internal_rotation_stop'& triallist$ppn ==ptest, ]
ex2 <- ex2[1, ]
#ex2$trialindex
exts2 <- read.csv(paste0(procd, 'zscal_pp_', ex2$ppn, '_trial_', ex2$trialindex, '_aligned.csv'))
exts2 <- dplyr::rename(exts2, rectus_abdominis = rectus_abdominus)
  exts2 <- subset(exts2, time_s_c > 1 & time_s_c <5)
  exts2 <- exts2[seq(1, nrow(exts2), by = 5), ]    #downsample by a 5th
  exts2$time_ms <- round(exts2$time_s_c*1000)
  fp <- pracma::findpeaks(exts2$speed)
  fp <- fp[which.max(fp[,1]),]
  startmov <- exts2$time_ms[fp[3]]
  endmov   <- exts2$time_ms[fp[4]]
  ps       <- exts2$time_ms[fp[2]]
#apply movement rule to centralize trials
#here all processing steps are performed, and data aligned
exampleenvelope <- ggplot(exts2, aes(x=time_ms)) +
  geom_path(aes(y=env), size =1) +
  theme_cowplot(12) + 
  ggtitle('Example trial internal rotation') +
  ylab('norm. env') +
  geom_vline(xintercept=startmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=endmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=startmov-500, color = 'black', size= 2, alpha=.3) +
  geom_smooth(aes(y=env), method='lm', linetype = 'dashed',color = 'black') +
  geom_vline(xintercept=ps, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept= 3000) +
  xlim(1000, 5000) +
  xlab('time (ms)')

examplemovement <- ggplot(exts2, aes(x=time_ms)) +
  geom_path(aes(y=round(speed*2500*100)), color = 'black') +
  theme_cowplot(12) +
  geom_vline(xintercept=startmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=endmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=startmov-500, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept=ps, color = 'black', size= 2, alpha=.3) +
  ylab("cm per second") +
  geom_vline(xintercept=3000, color = 'black', size= 1) +
  geom_vline(xintercept= 3000) +
  xlim(1000, 5000)+xlab('time (ms)')

exampleCOP <- ggplot(exts2, aes(x=time_ms)) +
  geom_path(aes(y=COPc)) +
  theme_cowplot(12) +
  geom_vline(xintercept=startmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=endmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=startmov-500, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept=ps, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept= 3000) +xlim(1000, 5000) +
  xlab('time (ms)')

exampleEMG <- ggplot(exts2, aes(x=time_ms)) +
  geom_path(aes(y=infraspinatus, color = 'infraspinatus')) +
  geom_path(aes(y=pectoralis_major, color = 'pectoralis major')) +
  geom_path(aes(y=rectus_abdominis, color = 'rectus abdominis')) +
  geom_path(aes(y=erector_spinae, color = 'erector spinae')) +
  geom_vline(xintercept=startmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=endmov, color = 'red', size= 2, alpha=.3) +
  geom_vline(xintercept=startmov-500, color = 'black', size= 2, alpha=.3) +
  geom_vline(xintercept=ps, color = 'black', size= 2, alpha=.3) +
  scale_color_manual(values = colors_mus) +
  labs(y="sEMG") +
  theme_cowplot(12) +
  theme(legend.position = c(0.8, 0.6)) +
  geom_vline(xintercept= 3000) +
  xlim(1000, 5000) +
  xlab('time (ms)')

Figure S3: Combined time series. One example trial and the associated signals are shown, for the internal rotation movement condition. At time = 0 s, the prompt is given to the participant to vocalize. We determine a detrending line using linear regression for the 1 to 5 s after the vocalization prompt. Note, at 3 s (3000 ms), there is a movement prompt. However, we determine our window where we assess peaks in signals at 500 ms before and after the movement onset/offset (using peakfinding function on the 2D speed time series of the wrist). In these trials, the analysis window is given in gray dashed bars, which is 500 ms after and before movement onset/offset.

This refers to Figure 1, panel C in the paper.

Show code setting up the data frame

overwrite = FALSE

if(overwrite)
{
tr_wd <- triallist
tr_wd$min_amp_c_around_move <- NA
tr_wd$min_amp_time_around_move <- NA
tr_wd$max_amp_c_around_move <- NA
tr_wd$max_amp_time_around_move  <- NA
tr_wd$max_pectoral <- NA
tr_wd$max_infra<- NA
tr_wd$max_rectus<- NA
tr_wd$max_erector<- NA
tr_wd$max_COPxc <- NA
tr_wd$max_COPyc <- NA
tr_wd$max_COPc  <- NA

#lets loop through the triallist and extract the files
tr_wd$uniquetrial <- paste0(tr_wd$trialindex, '_', triallist$ppn) 
  #remove duplicates so that we ignore repeated trials #NOTE in total we registered 25 repeats
howmany_repeated <- sum(duplicated(tr_wd$uniquetrial))
tr_wd <- tr_wd[!duplicated(tr_wd$uniquetrial),]

#make one large dataset
tsl <- data.frame()

for(ID in tr_wd$uniquetrial)
{
   #debugging
  #ID <- tr_wd$uniquetrial[27]
  tr_wd_s <- tr_wd[tr_wd$uniquetrial==ID,]
  
  ##load in row tr_wd
  pp <- tr_wd_s$ppn
  gender <- mtlist$sex[mtlist$ppn==pp]
  hand   <- mtlist$handedness[mtlist$ppn==pp]
  triali <- tr_wd_s$trialindex
  if(file.exists(paste0(procd, 'zscal_pp_', pp, '_trial_', triali, '_aligned.csv')))
  {
  ts <- read.csv(paste0(procd, 'zscal_pp_', pp, '_trial_', triali, '_aligned.csv'))
  ts$time_ms <- round(ts$time_s*1000)
  ts$time_ms <- ts$time_ms-min(ts$time_ms)
                                              #take the first 1 to 5 seconds
  ts <- ts[(ts$time_ms>1000) & (ts$time_ms < 5000), ]
  ts$env_z <- pracma::detrend(ts$env_z, tt = 'linear') # detrend the envelope

  #movement onset + - 500 milliseconds
  fp <- pracma::findpeaks(ts$speed)
  fp <- fp[which.max(fp[,1]),]
  startmov <- ts$time_ms[fp[3]]-500
  endmov   <- ts$time_ms[fp[4]]+500
  ts$movstart <- ts$time_ms-ts$time_ms[fp[3]] #get the time when the movement starts (which is the beginning of the rise to peak speed)
  
  subts <- ts[(ts$time_ms>startmov) & (ts$time_ms<endmov),] #subset
  subts$peaks <- NA
  indexspeed <- which(subts$time_ms==ts$time_ms[fp[2]])[1] #there are might be multiple adjacent values that have the peak speed
  subts$peaks[indexspeed] <- 'peak speed'
  indexend <- which(subts$time_ms==ts$time_ms[fp[4]])[1] #there are might be multiple adjacent values that have the peak speed
  subts$peaks[indexend] <- 'movement end'
  subts$movement_condition <- tr_wd_s$movement_condition
  subts$vocal_condition <- tr_wd_s$vocal_condition
  subts$weight_condition <- tr_wd_s$weight_condition
  subts$trial_number <-  tr_wd_s$trialindex
  
  #bind into big dataset
  tsl <- rbind(tsl, subts)

  #peakspeed
  tp <- subts$time_ms[which(subts$peaks=='peak speed')] #take everything before the peak speed so that we get movement onset
  
  #get the subts
  tr_wd$max_amp_c_around_move[tr_wd$uniquetrial==ID] <- max(subts$env_z[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  tr_wd$min_amp_c_around_move[tr_wd$uniquetrial==ID] <- min(subts$env_z[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  # also collect the times of the peaks
  envcc <- subts$env_z[subts$time_ms<tp]
  timecc <- subts$time_ms[subts$time_ms<tp]
  tr_wd$max_amp_time_around_move[tr_wd$uniquetrial==ID] <- timecc[which.max(envcc)]-tp # time pos peak
  tr_wd$min_amp_time_around_move[tr_wd$uniquetrial==ID] <- timecc[which.min(envcc)]-tp #time negative peak relative to peak speed
  # collect the peaks
  tr_wd$max_pectoral[tr_wd$uniquetrial==ID] <- max(subts$pectoralis_major[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  tr_wd$max_infra[tr_wd$uniquetrial==ID] <- max(subts$infraspinatus[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  tr_wd$max_rectus[tr_wd$uniquetrial==ID] <- max(subts$rectus_abdominis[subts$time_ms<tp], na.rm=TRUE) #before the peak speed
  tr_wd$max_erector[tr_wd$uniquetrial==ID] <- max(subts$erector_spinae[subts$time_ms<tp], na.rm=TRUE)#before the peak speed
  tr_wd$max_COPxc[tr_wd$uniquetrial==ID] <- max(subts$COPXc[subts$time_ms<tp], na.rm=TRUE)#before the peak speed
  tr_wd$max_COPyc[tr_wd$uniquetrial==ID] <- max(subts$COPYc[subts$time_ms<tp], na.rm=TRUE)#before the peak speed
  tr_wd$max_COPc[tr_wd$uniquetrial==ID] <- max(subts$COPc[subts$time_ms<tp], na.rm=TRUE)#before the peak speed
  }
}

#write a large dateset to 
#write.csv(tsl, )

#order levels
tr_wd$movement_condition <- recode(tr_wd$movement_condition, no_movement  = "no movement",
                                           extension_stop = "extension",flexion_stop  = "flexion",
                                           external_rotation_stop = "external rotation",
                                           internal_rotation_stop = "internal rotation")
tr_wd$movement_condition <- factor(tr_wd$movement_condition, 
                                   levels = c("no movement", "extension", "flexion", 
                                              "external rotation", "internal rotation"))

#maybe write to file?
write.csv(tr_wd, paste0(procdtot, 'fulldata.csv'))
write.csv(tsl, paste0(procdtot, 'time_series_all.csv'))
}
#read this in
if(!overwrite)
{
  tr_wd <- read.csv(paste0(procdtot, 'fulldata.csv'))
  tsl   <- read.csv(paste0(procdtot, 'time_series_all.csv'))
}

tr_wd$ppn <- ifelse(tr_wd$ppn%in%c(41, 42), 4, tr_wd$ppn)
tsl <- dplyr::rename(tsl, rectus_abdominis = rectus_abdominus)

Supporting results

Descriptive results

Fig. S4 shows an overview of the muscle activity patterns for each movement condition, split over weight conditions. Table S1 provides the numerical information. Table S2 provides normalized peak muscle activity by weight. All in all, these descriptive results pattern in sensible ways. Firstly, we can confirm that indeed the focal muscles powering internal (pectoralis major) or external (infraspinatus) rotation are peaking in activity during these movement conditions.

Show code producing the plot

library(magick)

##### This the raw code for the plot
    #we have however, done after editing to prettify and mark the plot
sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='vocalize',] #only real trials

sub$max_pectoral <- scale(sub$max_pectoral, center = FALSE)[,1]
sub$max_infra    <- scale(sub$max_infra, center = FALSE )[,1]
sub$max_rectus   <- scale(sub$max_rectus, center = FALSE)[,1]
sub$max_erector  <- scale(sub$max_erector, center = FALSE)[,1]
sub$max_COPc     <- scale(sub$max_COPc, center = FALSE)[,1]

#reorder values
sub$movement_condition <- factor(sub$movement_condition, levels=rev(c('no movement', 'internal rotation', 'external rotation', 'flexion', 'extension')))

peaks_pectoralis_bymovement <- ggplot(sub, aes(x = movement_condition, y=max_pectoral)) +
  geom_quasirandom(aes(color = weight_condition),alpha = 0.4, dodge.width = 0.8, cex = 1) +
  geom_boxplot(aes(group=interaction(movement_condition, weight_condition)),alpha = 0, size=0.5) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Pectoralis major (internal rotator)') +
  ylab('max sEMG activity \n (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = col_pectoralis)) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm')) +
  ylim(-0.1, 7) +
  geom_hline(yintercept=median(sub$max_pectoral[sub$movement_condition == 'no movement']), linetype='dashed') +
  coord_flip()

peaks_infraspinatus_bymovement <- ggplot(sub, aes(x = movement_condition,y=max_infra)) +
  geom_quasirandom(aes(color = weight_condition), alpha = 0.4, dodge.width = 0.8, cex = 1) +
  geom_boxplot(aes(group = interaction(movement_condition, weight_condition)), alpha = 0, size=0.5) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Infraspinatus (external rotator)') +
  ylab('max sEMG activity \n (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = col_infraspinatus)) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm')) +
  ylim(-0.1, 7) +
  geom_hline(yintercept = median(sub$max_infra[sub$movement_condition == 'no movement']), linetype='dashed') +
  coord_flip()

peaks_rectus_bymovement <- ggplot(sub, aes(x = movement_condition,y=max_rectus)) +
  geom_quasirandom(aes(color = weight_condition), alpha = 0.4, dodge.width = 0.8, cex = 1) +
  geom_boxplot(aes(group=interaction(movement_condition, weight_condition)), alpha = 0,  size=0.5) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Rectus abdominis (postural)') +
  ylab('max sEMG activity \n (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = col_rectus)) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm')) +
  ylim(-0.1, 7) +
  geom_hline(yintercept=median(sub$max_rectus[sub$movement_condition == 'no movement']), linetype='dashed') +
  coord_flip()

peaks_erector_bymovement <- ggplot(sub, aes(x = movement_condition, y=max_erector)) +
  geom_quasirandom(aes(color = weight_condition), alpha = 0.4, dodge.width = 0.8, cex = 1) +
  theme_cowplot(12) +
  geom_boxplot(aes(group=interaction(movement_condition, weight_condition)), alpha = 0, size=0.5) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Erector spinae (postural)')+
  ylab('max sEMG activity \n (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = col_erector)) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm') ) +
  ylim(-0.1, 7) +
  geom_hline(yintercept=median(sub$max_erector[sub$movement_condition == 'no movement']), linetype='dashed') +
  coord_flip()

# Center of Pressure
peaks_COP_bymovement <- ggplot(sub, aes(x = movement_condition,y= max_COPc)) +
  geom_quasirandom(aes(color = weight_condition), alpha = 0.4, dodge.width = 0.8, cex = 1) +
  theme_cowplot(12) +
  geom_boxplot(aes(group=interaction(movement_condition, weight_condition)), alpha = 0, size = 0.5) +
  theme(axis.text.x = element_text(angle = 0, size = 7)) +
  xlab("movement condition") +
  ggtitle('Change in center of pressure (COPc)') +
  ylab('center of pressure (normalized)') +
  theme(plot.title = element_text(hjust = 0.5, size = 10, color = 'grey')) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title=element_blank()) +
  theme(legend.position = c(0.8, 0.8), legend.key.size = unit(.2, 'cm')) +
  geom_hline(yintercept=median(sub$max_COPc[sub$movement_condition == 'no movement']),linetype='dashed') +
  coord_flip()

Figure S4: Manipulation check movement condition and peak muscle activity. The peak muscle activity, normalized for each muscle, is shown per movement and weight condition. A) highlights the fact that the pectoralis major, is - as to be expected - most active for the internal rotation movement condition. B) Also to be expected, the infraspinatus is most active during the external rotation. Note further, that both postural muscles seem more active during extension (and secondarily flexion).

This refers to Figure 1, panel D in the paper.

Show code producing the table

library(tidyr)
library(gmodels)
library(plyr)

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='vocalize',] #only real trials

sub$max_pectoral <- scale(sub$max_pectoral, center = FALSE)[,1]
sub$max_infra    <- scale(sub$max_infra, center = FALSE )[,1]
sub$max_rectus   <- scale(sub$max_rectus, center = FALSE)[,1]
sub$max_erector  <- scale(sub$max_erector, center = FALSE)[,1]
sub$max_COPc     <- scale(sub$max_COPc, center = FALSE)[,1]

subl <- gather(sub, key = muscle, value = peak_activity, c("max_pectoral", "max_infra", "max_rectus", "max_erector"))
subl$muscle <- revalue(subl$muscle, c("max_pectoral"="pectoralis major", "max_infra"="infraspinatus", "max_rectus"="rectus abdominis", "max_erector"="erector spinae"))

# #reorder values
# subl$movement_condition <- factor(subl$movement_condition, levels=rev(c('no movement', 'internal rotation', 'external rotation', 'flexion', 'extension')))

descriptivesbymusclesmovement <- subl %>%
  group_by(muscle, movement_condition) %>%
  dplyr::summarize(
    mean = round(ci(peak_activity)[1],2),
    sd = round(ci(peak_activity)[4],2),
    lowCI = round(ci(peak_activity)[2],2),
    hiCI = round(ci(peak_activity)[3], 2),)

#reorder to make it more comparable to previous Figure
muscle_order <- c("pectoralis major", "infraspinatus", "rectus abdominis", "erector spinae")
movement_order <- c("no movement", "internal rotation", "external rotation", "flexion", "extension")

descriptivesbymusclesmovement <- descriptivesbymusclesmovement %>%
  mutate(
    muscle = factor(muscle, levels = muscle_order),
    movement_condition = factor(movement_condition, levels = movement_order)
  )

#arrange table based on the defined order
descriptivesbymusclesmovement <- descriptivesbymusclesmovement %>%
  arrange(muscle, movement_condition)

#rename columns
descriptivesbymusclesmovement <- descriptivesbymusclesmovement %>% 
  dplyr::rename(`Muscle` = `muscle`,
                `Movement condition` = `movement_condition`,
                `Mean` = `mean`,
                `SD` = `sd`,
                `Low 95% CI` = `lowCI`,
                `High 95% CI` = `hiCI`)

Table S1: Normalized peak muscle activity for the different movement conditions. These are the numerical results associated with Fig. S4. The values indicate the peak EMG activity normalized for each muscle.

Muscle	Movement condition	Mean	SD	Low 95% CI	High 95% CI
pectoralis major	no movement	0.16	0.01	0.13	0.19
pectoralis major	internal rotation	1.09	0.10	0.90	1.28
pectoralis major	external rotation	0.56	0.03	0.50	0.62
pectoralis major	flexion	0.77	0.05	0.68	0.86
pectoralis major	extension	1.01	0.05	0.91	1.11
infraspinatus	no movement	0.04	0.00	0.03	0.05
infraspinatus	internal rotation	0.55	0.03	0.49	0.60
infraspinatus	external rotation	1.78	0.09	1.61	1.96
infraspinatus	flexion	0.32	0.02	0.29	0.35
infraspinatus	extension	0.35	0.02	0.31	0.40
rectus abdominis	no movement	0.81	0.05	0.71	0.91
rectus abdominis	internal rotation	0.89	0.04	0.81	0.96
rectus abdominis	external rotation	0.85	0.04	0.78	0.93
rectus abdominis	flexion	0.88	0.05	0.79	0.97
rectus abdominis	extension	0.94	0.04	0.86	1.03
erector spinae	no movement	0.37	0.02	0.33	0.41
erector spinae	internal rotation	0.61	0.03	0.56	0.67
erector spinae	external rotation	0.57	0.03	0.51	0.63
erector spinae	flexion	0.92	0.06	0.80	1.04
erector spinae	extension	1.32	0.08	1.16	1.48

This is referred to as Table S1 in the paper.

Further, postural muscles such as the rectus abdominis are especially active for extension movements (and secondarily flexion). Confirming their combined postural role, the muscle activity of the rectus abdominis and erector spinae are reliably correlated. Confirming their antagonistic role in rotating the humerus, the pectoralis and infraspinatus muscle activity are reliably correlated, indicative of their joint agonist/antagonistic control of posture and upper limb rotation, respectively. Lastly, from Table S2 it is clear that adding a wrist weight increases the peak muscle activity for all muscles, except the rectus abdominis (with no or barely overlapping 95 % confidence intervals in the weight vs. no weight condition).

Show code producing the table

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='vocalize',] #only real trials

sub$max_pectoral <- scale(sub$max_pectoral, center = FALSE)[,1]
sub$max_infra    <- scale(sub$max_infra, center = FALSE )[,1]
sub$max_rectus   <- scale(sub$max_rectus, center = FALSE)[,1]
sub$max_erector  <- scale(sub$max_erector, center = FALSE)[,1]
sub$max_COPc     <- scale(sub$max_COPc, center = FALSE)[,1]

subl <- gather(sub, key = muscle, value = peak_activity, c("max_pectoral", "max_infra", "max_rectus", "max_erector"))
subl$muscle <- revalue(subl$muscle, c("max_pectoral"="pectoralis major", "max_infra"="infraspinatus", "max_rectus"="rectus abdominis", "max_erector"="erector spinae"))
subl$weight_condition <- revalue(subl$weight_condition, c("no_weight"="no weight", "weight"="weight"))

descriptivesbymusclesweight <- subl %>%
  group_by(muscle, weight_condition) %>%
  dplyr::summarize(
    mean = round(ci(peak_activity)[1], 2),
    lowCI = round(ci(peak_activity)[2], 2),
    hiCI = round(ci(peak_activity)[3], 2),
    sd = round(ci(peak_activity)[4], 2))

#reorder to make it more comparable to previous Figure
muscle_order <- c("pectoralis major", "infraspinatus", "rectus abdominis", "erector spinae")
weight_order <- c("no weight", "weight")

descriptivesbymusclesweight <- descriptivesbymusclesweight %>%
  mutate(
    muscle = factor(muscle, levels = muscle_order),
    weight_condition = factor(weight_condition, levels = weight_order)
  )

#arrange table based on the defined order
descriptivesbymusclesweight <- descriptivesbymusclesweight %>%
  arrange(muscle, weight_condition)
            
#rename columns
descriptivesbymusclesweight <- descriptivesbymusclesweight %>% 
  dplyr::rename(`Muscle` = `muscle`,
                `Weight condition` = `weight_condition`,
                `Mean` = `mean`,
                `SD` = `sd`,
                `Low 95% CI` = `lowCI`,
                `High 95% CI` = `hiCI`)

#rearrange columns
descriptivesbymusclesweight <- descriptivesbymusclesweight %>%
  relocate(SD, .after = Mean)

Table S2: Normalized peak muscle activity per weight condition. The values indicate the peak EMG activity (normalized per muscle) per weight condition.

Muscle	Weight condition	Mean	SD	Low 95% CI	High 95% CI
pectoralis major	no weight	0.65	0.04	0.58	0.73
pectoralis major	weight	0.79	0.04	0.71	0.86
infraspinatus	no weight	0.54	0.04	0.47	0.61
infraspinatus	weight	0.69	0.05	0.59	0.79
rectus abdominis	no weight	0.87	0.03	0.81	0.92
rectus abdominis	weight	0.88	0.03	0.83	0.93
erector spinae	no weight	0.62	0.03	0.56	0.68
erector spinae	weight	0.90	0.04	0.83	0.98

This is referred to as Table S2 in the paper.

Exploratory visualization

Fig. S5 shows the amplitude envelope and EMG activity of a 1-second interval around the movement onset. This plot is informative about the rationale of the confirmatory analyses. Firstly, it can be seen that in general there is a positive peak during a movement onset in the amplitude envelope (relative to the trend line for that vocalization), that may or may not be preceded by a less extreme negative peak. Especially, for the internal rotation we observe such a negative peak. We will therefore assess whether particular movement conditions predict higher magnitude negative and positive peaks in the amplitude envelope (analysis 1). We further assess whether particular muscles predict the magnitude of positive or negative peaks (analysis 2). Finally, we will assess whether particular muscles are related to postural stability, to confirm the posture-mediating role of the muscles (analysis 3). In the current confirmatory analyses we focus on investigating the relationship between peaks, i.e. local minima and maxima in the different signals. This type of analysis is supported by a more straightforward power analysis than more complex continuous approaches (e.g. Generalized Additive Mixed Models, as underlying Fig. S5).

Show code producing the plot

#we will only look at vocalizations and we exclude practice trials, which have a number < 9
sub2 <- subset(tsl, vocal_condition != 'expire' & trial_number > 9)
sub2$movement_condition <- revalue(sub2$movement_condition, c(
  "extension_stop" = "extension",
  "internal_rotation_stop" = "internal r.",
  "flexion_stop" = "flexion",
  "external_rotation_stop" = "external r."))

#normalize
sub2$pectoralis_major  <- scale(sub2$pectoralis_major)
sub2$infraspinatus     <- scale(sub2$infraspinatus)
sub2$rectus_abdominis  <- scale(sub2$rectus_abdominis)
sub2$erector_spinae    <- scale(sub2$erector_spinae)

#what is the average movement start and end?
average_mov_start <- mean(sub2$time_ms[sub2$movstart==0])
#this we will use to create 

smoothenvelope <- ggplot(sub2, aes(x=movstart)) +
  geom_smooth(aes(y = env_z), method = 'gam', color = 'purple') +
  facet_grid(.~movement_condition)+geom_hline(yintercept = 0) +
  xlim(-500, 500) +
  theme_cowplot(12) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('amplitude envelope \n (normalized)') +
  xlab('time centered to movement onset')

smoothEMGactivity <- ggplot(sub2, aes(x = movstart)) +
  geom_smooth(aes(y = pectoralis_major, color = 'pectoralis major'), method = 'gam') +
  geom_smooth(aes(y = rectus_abdominis, color = 'rectus abdominis'), method = 'gam') +
  geom_smooth(aes(y = infraspinatus, color = 'infraspinatus'), method = 'gam') +
  geom_smooth(aes(y = erector_spinae, color = 'erector spinae'), method = 'gam') +
  xlim(-500, 500) +
  facet_grid(.~movement_condition) +
  theme_cowplot(12) +
  scale_color_manual(values = colors_mus) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('EMG activity \n (normalized)') +
  xlab('time centered to movement onset') +
  theme(legend.position = "none")    

smoothspeedCOP <- ggplot(sub2, aes(x = movstart)) +
  geom_smooth(aes(y = speed), color= 'orange', method = 'gam') +
  xlim(-500, 500) +
  facet_grid(.~movement_condition) +
  theme_cowplot(12) +
  scale_color_manual(values = colors_mus) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('speed \n (change in COP)') +
  xlab('time centered to movement onset')

smoothCOP <- ggplot(sub2, aes(x = movstart)) +
  geom_smooth(aes(y = COPc), method = 'gam', color = 'darkgrey') +
  xlim(-500, 500) +
  facet_grid(.~movement_condition) +
  theme_cowplot(12) +
  scale_color_manual(values = colors_mus) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
  ylab('change in center of pressure') +
  xlab('time centered to movement onset')

Figure S5: Smoothed conditional means for vocalization and muscle activity (with gam). Smoothed conditional means over time, where time is centered at 0 at the movement onset (as determined by the motion tracking of the wrist). The smoothed conditional means are generated by fitting non-linear smooths using generalized additive models in R-package ggplot2. It can be seen that especially for the extension movement there are clear anticipatory postural muscle activations (Aruin and Latash 1995), of the rectus abdominis before the movement onset, which is then followed by postural adjustments of the erector spinae. For the flexion condition, this activation pattern is reversed as one would expect given that the impulse vector should be directed in the opposite direction.

This is referred to as Figure S5 in the paper.

Main Confirmatory Analyses

As pre-registered, for the current confirmatory analyses we focus on the vocalization condition, ignoring the expiration baseline conditions as they are of secondary interest at this point of our inquiry.

Do different movement conditions have different effects on vocalization amplitude?

From inspecting the summarized trajectories, and the descriptive exploratory analyses above, we obtain that the internal rotation of the arm seems to start with a negative peak around the onset of the movement, which is followed by a positive peak. Furthermore, we observe that all other movements primarily have a positive peak. A straightforward test of whether the amplitude envelope has positive or negative peaks is to assess differences in peaks in the vocalization conditions per movement (and weight) condition.

Show code producing the plot

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='vocalize',] #only real trials

#reorder values & turn into factor
sub$movement_condition <- factor(sub$movement_condition, levels=c(
  "no movement",
  "internal rotation",
  "external rotation",
  "flexion",
  "extension"))
#rename no_weight in no weight
sub$weight_condition <- dplyr::recode(sub$weight_condition,
                               weight = "weight",
                               no_weight = "no weight")
#turn into factor
sub$weight_condition <- factor(sub$weight_condition, levels=c(
  "no weight",
  "weight"))


pospeakamplitudemovement <- ggplot(sub, aes(x = movement_condition, y = max_amp_c_around_move)) +
  geom_quasirandom(size = 2,alpha = 0.5, dodge.width=0.8, cex=1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("movement condition") +
  ggtitle("Positive peak amplitude") +
  ylab("amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title  =element_blank()) +
  geom_hline(yintercept = median(sub$max_amp_c_around_move[sub$movement_condition == "no movement"]), size=2, alpha=0.5, color = "blue") +
  ylim(-0.2, 2.3)

negpeakamplitudemovement <- ggplot(sub, aes(x = movement_condition, y = min_amp_c_around_move)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("movement condition") +
  ggtitle('Negative peak amplitude') +
  ylab("amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title = element_blank()) +
  geom_hline(yintercept = median(sub$min_amp_c_around_move[sub$movement_condition == "no movement"]), size = 2, alpha = 0.5, color = "blue") +
  ylim(-2.3, 0.2)

pospeakamplitudeweight <- ggplot(sub[sub$movement_condition != 'no movement', ], aes(x = weight_condition, y = max_amp_c_around_move)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("weight condition") +
  ggtitle("Positive peak amplitude") +
  ylab("amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title=element_blank()) +
  ylim(-0.2, 2.3)

negpeakamplitudeweight <- ggplot(sub[sub$movement_condition != 'no movement', ], aes(x = weight_condition,y = min_amp_c_around_move)) +
  geom_quasirandom(size = 2, alpha = 0.5, dodge.width = 0.8, cex = 1, color = "#476378ff") +
  geom_boxplot(alpha = 0, linewidth = 0.75) +
  theme_cowplot(12) +
  theme(axis.text.x = element_text(angle = 90, size = 10)) +
  xlab("weight condition") +
  ggtitle("Negative peak amplitude") +
  ylab("amplitude") +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.title=element_blank()) +
  ylim(-2.3, 0.2)

Figure S6: Effects of movement condition of positive and negative peaks in the voice amplitude. The upper part shows the positive peaks in the amplitude envelope during the different movement and weight conditions. The lower part shows the negative peaks (hence the negative values; note, in the modeling we will absolutize these values). It is clear that relative to vocalization of no movement, there are especially positive, but also more negative peaks in the amplitude envelope for the different movement conditions.

This is referred to as Figure S6 in the paper.

Show code running models on positive amplitude peaks

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='vocalize',] #only real trials
sub$movement_condition <- factor(sub$movement_condition, levels=c('no movement', 'internal rotation', 'external rotation', 'flexion', 'extension'))
sub$weight_condition <- dplyr::recode(sub$weight_condition,
                               weight = "weight",
                               no_weight = "no weight")
sub$weight_condition <- factor(sub$weight_condition, levels=c("no weight", "weight"))

#model 1
pospeaks_nullmodel <- lmer(max_amp_c_around_move ~ 1+(1|ppn), data = sub)
pospeaks_model   <- lmer(max_amp_c_around_move ~ weight_condition + movement_condition + (1|ppn), data = sub)
pospeaks_comparison <- anova(pospeaks_nullmodel, pospeaks_model)

#table matrix
  #Estimate, 

modelsummarypospeaksweightmovement <- summary(pospeaks_model)$coefficients
colnames(modelsummarypospeaksweightmovement) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarypospeaksweightmovement) <- c(
  'Intercept (no movement/no weight)',
  'vs. weight',
  'vs. internal rotation',
  'vs. external rotation',
  'vs. flexion',
  'vs. extension')

#posthoc (not necessary at this point)
#posthoc <- emmeans(modelmax, list(pairwise ~ movement_condition), adjust = "tukey")
#pcomparisons <- posthoc$`pairwise differences of movement_condition`

  #########################################also make a comparison with only movement conditions for the weight condition (CHANGECONF)
subsub <- subset(sub, movement_condition != 'no movement')
#model 1
pospeaks_nullmodel_onlymovement <- lmer(max_amp_c_around_move ~ 1 + (1 | ppn), data = subsub)
pospeaks_model_onlymovement   <- lmer(max_amp_c_around_move ~ weight_condition + (1 | ppn), data = subsub)
pospeaks_comparison_onlymovement <- anova(pospeaks_nullmodel_onlymovement, pospeaks_model_onlymovement)

#table matrix
  #Estimate, 

modelsummarypospeaksweightmovement_onlymovement <- summary(pospeaks_model_onlymovement)$coefficients
colnames(modelsummarypospeaksweightmovement_onlymovement) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarypospeaksweightmovement_onlymovement) <- c('Intercept (no weight)', 'vs. weight')

We first modeled with a mixed linear regression the variation in positive peaks in the amplitude envelope (using R-package lme4), with participant as random intercept (for more complex random slope models did not converge). A model with weight and movement conditions explained more variance than a base model predicting the overall mean, Change in Chi-Squared (5.00) = 36.06, p < .001). The model coefficients are given in Table S3. It can be seen that there is a positive but not statistically reliable effect of wrist weight in this sample. Further, all movements (extension, flexion, internal rotation, external rotation) lead to statistically reliable increases in positive peaks in the amplitude envelope relative to the no movement condition (with flexion and external rotation leading to more extreme effects).

Table S3: Effects of weight and movement condition on magnitude positive peaks in amplitude envelope.

	Estimate	SE	df	t-value	p-value
Intercept (no movement/no weight)	0.233	0.037	42.576	6.361	0.000
vs. weight	0.058	0.021	603.790	2.825	0.005
vs. internal rotation	0.104	0.033	602.729	3.204	0.001
vs. external rotation	0.116	0.032	602.740	3.586	0.000
vs. flexion	0.164	0.033	602.778	5.057	0.000
vs. extension	0.130	0.032	602.718	4.005	0.000

This is referred to as Table S3 in the paper.

Note that in the previous model we assess the role of weight vs. no weight. But during the no movement condition the weight is not affecting the participant. Thus, a better estimate of the effect of wrist weight is to assess weight vs. no weight conditions for only the movement conditions. Below it shows indeed that this comparison further increases the effect strength as compared to a model with the no movement condition included.

Table S4: Effects of weight on magnitude of positive peaks in amplitude envelope (no movement condition filtered out).

	Estimate	SE	df	t-value	p-value
Intercept (no weight)	0.361	0.035	20.170	10.296	0.00
vs. weight	0.062	0.024	484.657	2.570	0.01

This is referred to as Table S4 in the paper.

Show code running models on negative amplitude peaks

# model 2
sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='vocalize',] #only real trials
sub$movement_condition <- factor(sub$movement_condition, levels=c('no movement', 'internal rotation', 'external rotation', 'flexion', 'extension'))
sub$weight_condition <- dplyr::recode(sub$weight_condition,
                               weight = "weight",
                               no_weight = "no weight")
sub$weight_condition <- factor(sub$weight_condition, levels=c("no weight", "weight"))

sub$min_amp_c_around_move_abs <- abs(sub$min_amp_c_around_move)

negpeaks_nullmodel <- lmer(min_amp_c_around_move_abs ~  1  + (1|ppn), data = sub)
negpeaks_model <- lmer(min_amp_c_around_move_abs ~  weight_condition + movement_condition + (1|ppn), data = sub)
negpeaks_comparison <- anova(negpeaks_nullmodel, negpeaks_model)


#table matrix
modelsummarynegpeaksweightmovement <- summary(negpeaks_model)$coefficients
colnames(modelsummarynegpeaksweightmovement) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarynegpeaksweightmovement) <- c(
'Intercept (no movement/no weight)',
'vs. weight',
'vs. internal rotation',
'vs. external rotation',
'vs. flexion',
'vs. extension')

#also make a comparison with only movement conditions for the weight condition (CHANGECONF)
subsub <- subset(sub, movement_condition != 'no movement')

negpeaks_nullmodel_onlymovement <- lmer(min_amp_c_around_move ~ 1 + (1 | ppn), data = subsub)
negpeaks_model_onlymovement   <- lmer(min_amp_c_around_move ~ weight_condition + (1 | ppn), data = subsub)
negpeaks_comparison_onlymovement <- anova(negpeaks_nullmodel_onlymovement, negpeaks_model_onlymovement)

modelsummarynegpeaksweightmovement_onlymovement <- summary(negpeaks_model_onlymovement)$coefficients
colnames(modelsummarynegpeaksweightmovement_onlymovement) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarynegpeaksweightmovement_onlymovement) <- c('Intercept (no weight)', 'vs. weight')

Secondly, we modeled in a similar way the negative peaks in the amplitude envelope, and found that a model with weight and movement conditions explained more variance than a base model predicting the overall mean, Change in Chi-Squared (5.00) = 48.54, p < .001). Model coefficients are shown in Table S5. We find that some movement conditions (flexion, internal rotation, and especially external rotation) had higher magnitude negative peaks relative to no movement. No reliable effect of weight condition was found, nor did the extension movement lead to negative peaks relative to the no movement condition. Note, in our models, we have absolutized the values of negative peaks, such that positive effects of some condition means higher magnitude negative peaks (i.e., more negative peaks).

Table S5: Effects of weight and movement condition on magnitude negative peaks in amplitude envelope.

	Estimate	SE	df	t-value	p-value
Intercept (no movement/no weight)	0.185	0.024	30.928	7.747	0.000
vs. weight	0.020	0.011	603.514	1.784	0.075
vs. internal rotation	0.110	0.018	602.850	6.102	0.000
vs. external rotation	0.094	0.018	602.857	5.238	0.000
vs. flexion	0.094	0.018	602.879	5.243	0.000
vs. extension	0.080	0.018	602.843	4.478	0.000

This is referred to as Table S5 in the paper.

Table S6: Effects of weight on magnitude of negative peaks in amplitude envelope (no movement condition filtered out).

	Estimate	SE	df	t-value	p-value
Intercept (no weight)	-0.274	0.024	18.500	-11.611	0.000
vs. weight	-0.029	0.013	484.464	-2.264	0.024

This is referred to as Table S6 in the paper.

Is different muscle activity differently affecting amplitude of vocalization?

Show code calculating VIF

sub <- tr_wd[tr_wd$vocal_condition == 'vocalize' & tr_wd$trialindex>9,]
#check VIF
library(usdm)
df = cbind.data.frame(sub$max_COPc, sub$max_erector, sub$max_infra, sub$max_pectoral, sub$max_rectus) # Data Frame
VIFFIES <- vif(df)

Since each movement and weight condition is designed to recruit different muscle activations, we can also directly relate muscle activity peaks with the positive and negative peaks in the amplitude envelope. We use a similar linear mixed regression approach to model variance in vocal amplitude peaks with peaks in muscle activity for the different muscles measured. Since there are correlations between these muscle activities (see Tables S10, S11, S12, S13), we first assessed the VIF’s between the muscle activity peaks, which yielded a maximum VIF value of 1.29. Since this is considered a low value (VIF > 5 is generally considered problematic), we can combine the different muscle activity measurements in one model to predict amplitude envelope peaks. Figure S7 shows the graphical results of these relationships.

Show code producing the plot

sub <- tr_wd[tr_wd$vocal_condition == 'vocalize' & tr_wd$trialindex>9,]

#normalize
sub$max_pectoral  <- scale(sub$max_pectoral, center = FALSE)
sub$max_infra     <- scale(sub$max_infra, center = FALSE)
sub$max_rectus    <- scale(sub$max_rectus, center = FALSE)
sub$max_erector   <- scale(sub$max_erector, center = FALSE)


pospeakamplitude_peakpectoralis <- ggplot(sub, aes(y = abs(max_amp_c_around_move))) +
  geom_point(aes(x = max_pectoral, color = 'pectoralis major'), size =0.2) +
  geom_smooth(aes(x = max_pectoral, color = 'pectoralis major'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('positive peak \n amplitude envelope') +
  xlab('peak sEMG \n pectoralis major') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7)

pospeakamplitude_peakrectus <- ggplot(sub, aes(y = abs(max_amp_c_around_move))) +
  geom_point(aes(x = max_rectus, color = 'rectus abdominis'), size = 0.2) +
  geom_smooth(aes(x = max_rectus, color = 'rectus abdominis'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('positive peak \n amplitude envelope') +
  xlab('peak sEMG \n rectus abdominis') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7)

pospeakamplitude_peakerector <- ggplot(sub, aes(y=abs(max_amp_c_around_move))) +  
  geom_point(aes(x = max_erector, color = 'erector spinae'), size =0.2) +
  geom_smooth(aes(x = max_erector, color = 'erector spinae'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('positive peak \n amplitude envelope') +
  xlab('peak sEMG \n erector spinae') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7)

pospeakamplitude_peakinfraspinatus <- ggplot(sub, aes(y = abs(max_amp_c_around_move))) +
  geom_point(aes(x = max_infra, color = 'infraspinatus'), size =0.2) +
  geom_smooth(aes(x = max_infra, color = 'infraspinatus'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('positive peak \n amplitude envelope') +
  xlab('peak sEMG \n infraspinatus') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7)


negpeakamplitude_peakpectoralis <- ggplot(sub, aes(y = min_amp_c_around_move)) +  
  geom_point(aes(x = max_pectoral, color = 'pectoralis major'), size =0.2) +
  geom_smooth(aes(x = max_pectoral, color = 'pectoralis major'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('negative peak \n amplitude envelope') +
  xlab('peak sEMG \n pectoralis major') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7)

negpeakamplitude_peakrectus <- ggplot(sub, aes(y = min_amp_c_around_move)) +  
  geom_point(aes(x = max_rectus, color = 'rectus abdominis'), size =0.2) +
  geom_smooth(aes(x = max_rectus, color = 'rectus abdominis'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position="none") +
  ylab('negative peak \n amplitude envelope') +
  xlab('peak sEMG \n rectus abdominis') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7)

negpeakamplitude_peakerector <- ggplot(sub, aes(y = min_amp_c_around_move)) +
  geom_point(aes(x = max_erector, color = 'erector spinae'), size = 0.2) +
  geom_smooth(aes(x = max_erector, color = 'erector spinae'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('negative peak \n amplitude envelope') +
  xlab('peak sEMG \n erector spinae') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7)

negpeakamplitude_peakeinfraspinatus <- ggplot(sub, aes(y = min_amp_c_around_move)) +
  geom_point(aes(x = max_infra, color = 'infraspinatus'), size =0.2) +
  geom_smooth(aes(x = max_infra, color = 'infraspinatus'), method = 'lm') +
  scale_color_manual(values = colors_mus) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab('negative peak \n amplitude envelope') +
  xlab('peak sEMG \n infraspinatus') +
  theme_cowplot() +
  theme(legend.position = "none") +
  xlim(0, 7)

This refers to Figure 2, panel C (positive) and panel D (negative) of the paper.

Figure S7: Relations between peak muscle activity and positive peaks in the amplitude envelope

Show code modeling positive amplitude peaks by muscle

sub <- tr_wd[ tr_wd$trialindex>9 & tr_wd$vocal_condition=='vocalize',] #only real trials

#model 1
pospeaks_bymuscle_nullmodel <- lmer(max_amp_c_around_move ~ 1 + (1|ppn), data = sub)
pospeaks_bymuscle <- lmer(max_amp_c_around_move ~ max_erector + max_infra + max_pectoral + max_rectus + (1 | ppn), data = sub)
pospeaks_bymuscle_comparison <- anova(pospeaks_bymuscle_nullmodel, pospeaks_bymuscle)

#get model summary
modelsummarypospeaksbymuscle <- summary(pospeaks_bymuscle)$coefficients
colnames(modelsummarypospeaksbymuscle) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarypospeaksbymuscle) <- c('Intercept', 'erector spinae', 'infraspinatus', 'pectoralis major', 'rectus abdominis')

A model with participant as random intercept (a more complex random slope model did not converge) the different peak muscle activities explained more variance than a base model predicting the overall mean, Change in Chi-Squared = 59.21 (4.00), p < .001). The model coefficients are given in Table S7. Further, Table S7 shows that peak EMG activity in all the muscles (except for the rectus abdominis) leads to statistically reliable increases in positive peaks in the amplitude envelope (with stronger effects).

Table S7: Linear mixed regression model assessing the relation of peak muscle activity with the positive peak in the amplitude envelope.

	Estimate	SE	df	t-value	p-value
Intercept	0.240	0.037	54.499	6.429	0.000
erector spinae	0.016	0.004	617.693	4.207	0.000
infraspinatus	0.006	0.002	612.279	2.832	0.005
pectoralis major	0.013	0.003	617.423	4.819	0.000
rectus abdominis	-0.003	0.007	612.589	-0.371	0.711

This is referred to as Table S7 in the paper.

We similarly modeled the variance of the magnitude in negative peaks in the amplitude envelope with muscle activity peaks (see right side of Figure S7 for graphical results), and we similarly observed that such a model performed better than a base model predicting the overall mean, Change in Chi-Squared = 16.48 (4.00), p = 0.002. The model coefficients are given in Table S8. Further, as shown in Table S8 for the infraspinatus and the pectoralis major, increases in peak EMG activity lead to statistically reliable increases in the magnitude of negative peaks in the amplitude envelope (with stronger effects), while the effects of erector spinae and the rectus abdominis do not reach statistical reliability.

Table S8: Linear mixed regression model assessing the relation of peak muscle activity with the negative peak in the amplitude envelope. The negative peaks in the amplitude envelope are absolutized.

	Estimate	SE	df	t-value	p-value
Intercept	0.227	0.026	36.153	8.830	0.000
erector spinae	0.004	0.002	613.122	1.943	0.053
infraspinatus	0.003	0.001	608.863	2.328	0.020
pectoralis major	0.003	0.002	612.841	2.152	0.032
rectus abdominis	0.001	0.004	619.957	0.359	0.719

This is referred to as Table S8 in the paper.

Is a particular muscle related to postural stability?

Show code modeling change in posture by muscle

sub <- tr_wd[tr_wd$vocal_condition == 'vocalize' & tr_wd$trialindex>9, ]
sub$max_COPc <- scale(sub$max_COPc)

maxCOPc_bymuscle_nullmodel <- lmer(max_COPc ~ 1 + (1 | ppn), data = sub)
maxCOPc_bymuscle <- lmer(max_COPc ~ max_erector + max_infra + max_pectoral + max_rectus + (1 | ppn), data = sub)
maxCOPc_bymuscle_comparison <- anova(maxCOPc_bymuscle_nullmodel, maxCOPc_bymuscle)

#get model summary and prepare table
modelsummarymaxCOPcbymuscle <- summary(maxCOPc_bymuscle)$coefficients
colnames(modelsummarymaxCOPcbymuscle) <- c('Estimate', 'SE', 'df', 't-value', 'p-value')
rownames(modelsummarymaxCOPcbymuscle) <- c('Intercept', 'erector spinae', 'infraspinatus', 'pectoralis major', 'rectus abdominis')

Finally, we will assess which muscle activity can be related to changes in the center of pressure, which would directly confirm that gesture-speech biomechanics relate to postural stability (Pouw, Harrison, and Dixon 2020). Figure S8 shows the graphical results. We similarly performed a linear mixed regression (with participant as random intercept) with a model containing peak EMG activity for each muscle which was regressed on the peak in change in the center of pressure (COPc). We obtained that a base model predicting the over all mean of COPc was outperformed relative to said model, Change in Chi-Squared = 180.01 (4.00), p < .001). Table S9 provides the coefficient information. We find that only the postural muscles (rectus abdominis, erector spinae) indeed reliably predict the magnitude of changes in the center of pressure, while the pectoralis major and the infraspinatus do not reliably relate to the changes in ground reaction forces.

Show code producing the plot

sub <- tr_wd[ tr_wd$trialindex>9,] #only real trials

sub$max_pectoral <- scale(sub$max_pectoral)
sub$max_infra    <- scale(sub$max_infra )
sub$max_rectus   <- scale(sub$max_rectus)
sub$max_erector  <- scale(sub$max_erector)
sub <- subset(sub, movement_condition != 'no movement' & vocal_condition == 'vocalize')

#correlation plot, to confirm role of postural muscle activations
COPcbymuscle_infraspinatus <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_infra, color = 'infraspinatus'), size = 0.2) +
  geom_smooth(aes(y = max_infra, color = 'infraspinatus'), method='lm') +
  ggtitle('Focal muscles: \n infraspinatus') +
  scale_color_manual(values = colors_mus) +
  ylab('max sEMG activity \n (normalized)') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 5, color = col_infraspinatus)) +
  theme(legend.position = "none") +
  theme_cowplot() +
  theme(legend.position = "none") #+facet_grid(.~movement_condition)

COPcbymuscle_pectoralis <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_pectoral, color = 'pectoralis major'), size=0.2) +
  geom_smooth(aes(y = max_pectoral, color = 'pectoralis major'), method = 'lm') +
  ggtitle('Focal muscles: \n pectoralis major') +
  scale_color_manual(values = colors_mus) +
  ylab('max sEMG activity \n (normalized)') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 5, color = col_pectoralis)) +
  theme(legend.position = "none") +
  theme_cowplot() +
  theme(legend.position = "none")  #+facet_grid(.~movement_condition)

COPcbymuscle_rectus <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_rectus, color = 'rectus abdominis'), size = 0.2) +
  geom_smooth(aes(y = max_rectus, color = 'rectus abdominis'), method = 'lm') +
  #ggtitle('Postural muscles') +
  scale_color_manual(values = colors_mus) +
  ylab('max sEMG activity \n (normalized)') +
  ggtitle('Postural muscles: \n rectus abdominis') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 5, color = col_rectus)) +
  theme(legend.position = "none") +
  theme_cowplot() +
  theme(legend.position = "none")  #+facet_grid(.~movement_condition)


COPcbymuscle_erector <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_erector, color = 'erector spinae'), size = 0.2) +
  geom_smooth(aes(y = max_erector, color = 'erector spinae'), method = 'lm') +
  #ggtitle('Postural muscles') +
  scale_color_manual(values = colors_mus) +
  ylab('max sEMG activity \n (normalized)') +
  ggtitle('Postural muscles: \n erector spinae') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 5, color = col_erector)) +
  theme(legend.position = "none") +
  theme_cowplot() +
  theme(legend.position = "none")  #+facet_grid(.~movement_condition)

The bottom-right plot (green) can be found in Figure 2, panel E.

Figure S8: Confirmation of postural muscle activation during changes in center of mass.

Table S9: Linear mixed regression model for predicting peak change in center of pressure based on muscle activity.

	Estimate	SE	df	t-value	p-value
Intercept	-0.875	0.122	55.418	-7.169	0.000
erector spinae	0.118	0.012	617.831	9.752	0.000
infraspinatus	0.028	0.007	612.430	4.168	0.000
pectoralis major	0.063	0.009	617.566	7.065	0.000
rectus abdominis	0.025	0.023	611.990	1.109	0.268

This is referred to as Table S9 in the paper.

Other exploratory analyses

Along with the confirmatory analyses reported above that relate to questions posed in the pre-registration, we report on additional, exploratory analyses. Firstly we assess for possible context-dependent interrelationships of posture, muscle activity, and voice, by producing correlation matrices per movement condition (see Tables S10, S11, S12, S13) which are represented in our synergies plot in Fig. 2 (panel F) in the main article. We see that producing different movements can recruit variable interrelationships between muscle, posture and voice. Striking to pick out is for example that we find a direct relationship between the degree of center of pressure change and the magnitude of positive amplitude peaks in the voice.

Show code producing the correlation matrices and plots

library(corx)
#correlation between positive and negative peaks in amplitude and change in center of pressure
correlation_pospeakCOPc_infraspinatus <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = max_amp_c_around_move, color = 'infraspinatus'), size = 0.2) +
  geom_smooth(aes(y = max_amp_c_around_move), method = 'lm') +
  ggtitle('Focal muscles: \n infraspinatus') +
  scale_color_manual(values = colors_mus) +
  ylab('positive amplitude peaks') +
  xlab('max. change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 20, color = col_infraspinatus),
        legend.position = "none") #+facet_grid(.~movement_condition)

correlation_negpeakCOPc_infraspinatus <- ggplot(sub, aes(x = max_COPc)) +
  geom_point(aes(y = min_amp_c_around_move, color = 'infraspinatus'), size = 0.2) +
  geom_smooth(aes(y = min_amp_c_around_move), method = 'lm') +
  ggtitle('Focal muscles: \n infraspinatus') +
  scale_color_manual(values = colors_mus) +
  ylab('negative amplitude peaks') +
  xlab('max change in \n center of pressure') +
  theme_cowplot(12) +
  theme(plot.title = element_text(hjust = 0.5, size = 20, color = col_infraspinatus),
        legend.position = "none") #+facet_grid(.~movement_condition)

#internal rotation
condition_internalrotation <- sub$movement_condition=='internal rotation'
data_condition_internalrotation <- cbind.data.frame(sub$max_pectoral[condition_internalrotation], sub$max_infra[condition_internalrotation], sub$max_rectus[condition_internalrotation], sub$max_erector[condition_internalrotation], sub$max_COPc[condition_internalrotation], sub$max_amp_c_around_move[condition_internalrotation])

#rename columns
data_condition_internalrotation_rn <- data_condition_internalrotation %>% 
  dplyr::rename(
    `max. EMG pectoralis major` = `sub$max_pectoral[condition_internalrotation]`,
    `max. EMG infraspinatus` = `sub$max_infra[condition_internalrotation]`,
    `max. EMG rectus abdominis` = `sub$max_rectus[condition_internalrotation]`,
    `max. EMG erector spinae` = `sub$max_erector[condition_internalrotation]`,
    `max. change in center of pressure` = `sub$max_COPc[condition_internalrotation]`,
    `max. positive amplitude vocalization` = `sub$max_amp_c_around_move[condition_internalrotation]`
    )
    
corr_internalrotation <- data_condition_internalrotation_rn %>%
  corx(triangle = "lower",
       stars    = c(0.05, 0.01, 0.001),
       note = "Internal rotation") 
#knitr::kable(corr_internalrotation$apa)


#external rotation
condition_externalrotation <- sub$movement_condition=='external rotation'
data_condition_externalrotation <- cbind.data.frame(sub$max_pectoral[condition_externalrotation], sub$max_infra[condition_externalrotation], sub$max_rectus[condition_externalrotation], sub$max_erector[condition_externalrotation], sub$max_COPc[condition_externalrotation], sub$max_amp_c_around_move[condition_externalrotation])

#rename columns
data_condition_externalrotation_rn <- data_condition_externalrotation %>% 
  dplyr::rename(
    `max. EMG pectoralis major` = `sub$max_pectoral[condition_externalrotation]`,
    `max. EMG infraspinatus` = `sub$max_infra[condition_externalrotation]`,
    `max. EMG rectus abdominis` = `sub$max_rectus[condition_externalrotation]`,
    `max. EMG erector spinae` = `sub$max_erector[condition_externalrotation]`,
    `max. change in center of pressure` = `sub$max_COPc[condition_externalrotation]`,
    `max. positive amplitude vocalization` = `sub$max_amp_c_around_move[condition_externalrotation]`
    )

corr_externalrotation <- data_condition_externalrotation_rn %>%
  corx(triangle = "lower",
       stars    = c(0.05, 0.01, 0.001),
       note = "External rotation") 
#knitr::kable(corr_externalrotation$apa)


#extension
condition_extension <- sub$movement_condition=='extension'
data_condition_extension <- cbind.data.frame(sub$max_pectoral[condition_extension], sub$max_infra[condition_extension], sub$max_rectus[condition_extension], sub$max_erector[condition_extension], sub$max_COPc[condition_extension], sub$max_amp_c_around_move[condition_extension])
#cor.test(sub$max_rectus[condition_extension], sub$max_COPc[condition_extension])

#rename columns
data_condition_extension_rn <- data_condition_extension %>% 
  dplyr::rename(
    `max. EMG pectoralis major` = `sub$max_pectoral[condition_extension]`,
    `max. EMG infraspinatus` = `sub$max_infra[condition_extension]`,
    `max. EMG rectus abdominis` = `sub$max_rectus[condition_extension]`,
    `max. EMG erector spinae` = `sub$max_erector[condition_extension]`,
    `max. change in center of pressure` = `sub$max_COPc[condition_extension]`,
    `max. positive amplitude vocalization` = `sub$max_amp_c_around_move[condition_extension]`
    )

corr_extension <- data_condition_extension_rn %>%
  corx(triangle = "lower",
       stars    = c(0.05, 0.01, 0.001),
       note = "Extension")
#knitr::kable(corr_extension$apa)


#flexion
condition_flexion <- sub$movement_condition=='flexion'
data_condition_flexion <- cbind.data.frame(sub$max_pectoral[condition_flexion], sub$max_infra[condition_flexion], sub$max_rectus[condition_flexion], sub$max_erector[condition_flexion], sub$max_COPc[condition_flexion], sub$max_amp_c_around_move[condition_flexion])

#rename columns
data_condition_flexion_rn <- data_condition_flexion %>% 
  dplyr::rename(
    `max. EMG pectoralis major` = `sub$max_pectoral[condition_flexion]`,
    `max. EMG infraspinatus` = `sub$max_infra[condition_flexion]`,
    `max. EMG rectus abdominis` = `sub$max_rectus[condition_flexion]`,
    `max. EMG erector spinae` = `sub$max_erector[condition_flexion]`,
    `max. change in center of pressure` = `sub$max_COPc[condition_flexion]`,
    `max. positive amplitude vocalization` = `sub$max_amp_c_around_move[condition_flexion]`
    )

corr_flexion <- data_condition_flexion_rn %>%
  corx(triangle = "lower",
       stars    = c(0.05, 0.01, 0.001),
       note = "Flexion") 
#knitr::kable(corr_flexion$apa)


#very stable correlation of posture and internal rotation
correlation_posture_internalrotation <- sub %>% 
  dplyr::filter(movement_condition == "internal rotation") %>% 
  ggplot(aes(x = max_COPc, y = max_amp_c_around_move)) +
  geom_point(size = 3) +
  geom_smooth(method = 'lm', color = 'red', size = 2) +
  theme_cowplot() +
  xlab('max change in \n center of pressure') +
  ylab('positive\namplitude peaks')
  
# ggplot(sub[sub$movement_condition=='internal rotation', ], aes(x = max_COPc, y = max_amp_c_around_move))+geom_point(size = 3)+geom_smooth(method='lm', color = 'red', size = 2)+theme_cowplot()+xlab('max change in \n center of pressure')+ylab('positive\namplitude peaks')

Table S10 shows the correlation matrix for the maximum peaks in EMG (row and columns 1-4), change in center of pressure (5) and maximum vocal amplitude (6). For this table, the data only include the internal rotation condition.

Tables S10, S11, S12, S13 refer to to the synergies plot in Fig. 2 (F) in the paper.

Table S10: Correlation matrix for the internal rotation condition.

	1	2	3	4	5
max. EMG pectoralis major
max. EMG infraspinatus	.34***
max. EMG rectus abdominis	.00	-.16
max. EMG erector spinae	.17	.09	.16
max. change in center of pressure	.30***	.23*	-.07	.21*
max. positive amplitude vocalization	.42***	.22*	-.02	.10	.51***

Table S11 shows the correlation matrix for the maximum peaks in EMG (row and columns 1-4), change in center of pressure (5) and maximum vocal amplitude (6) For this table, the data only include the external rotation condition.

Table S11: Correlation matrix for the external rotation condition.

	1	2	3	4	5
max. EMG pectoralis major
max. EMG infraspinatus	.04
max. EMG rectus abdominis	-.14	.06
max. EMG erector spinae	.20*	.07	.15
max. change in center of pressure	.12	.12	-.09	.27**
max. positive amplitude vocalization	.27**	.18*	-.05	.25**	.04

Table S12 shows the correlation matrix for the maximum peaks in EMG (row and columns 1-4), change in center of pressure (5) and maximum vocal amplitude (6) For this table, the data only include the extension condition.

Table S12: Correlation matrix for the extension condition.

	1	2	3	4	5
max. EMG pectoralis major
max. EMG infraspinatus	.18*
max. EMG rectus abdominis	.06	-.02
max. EMG erector spinae	.20*	.05	.15
max. change in center of pressure	.21*	.32***	.18*	.17
max. positive amplitude vocalization	-.02	.29***	-.01	.26**	.17

Table S13 shows the correlation matrix for the maximum peaks in EMG (row and columns 1-4), change in center of pressure (5) and maximum vocal amplitude (6) For this table, the data only include the flexion condition.

Table S13: Correlation matrix for the flexion condition.

	1	2	3	4	5
max. EMG pectoralis major
max. EMG infraspinatus	.18*
max. EMG rectus abdominis	-.19*	-.02
max. EMG erector spinae	.03	.26**	.05
max. change in center of pressure	.13	.34***	.02	.04
max. positive amplitude vocalization	.09	.23**	.00	.07	.16

Figure S9: Correlation of posture and amplitude for internal rotation.

This refers to Figure 2, panel G in the paper.

Figure S9 visualizes the positive correlation between posture (max. change in center of pressure) and vocalization amplitude (positive peaks in amplitude) found for the internal rotation condition in Table S10.

We finally explored the timing relationship of negative and positive peaks in the amplitude of the voice relative to the peak in speed of the hand. The density plot in Fig. S10 visualizes the temporal relationship between peak speed in movement and negative (left) and positive (right) peaks in the amplitude envelope, relative to the moment the moving hand reaches a peak in speed during the trial (i.e., t = 0, at peak speed of the wrist). Our findings corroborate our general observation that positive (rather than negative) peaks seem to be more robustly driven by gesture kinetics. In the timing relations obtained, the negative peaks that we observe are more variably distributed around peak speed, while for the positive peaks we see clear occurrences just around the movement onset (just before the peak speed at t = 0). This result, next to our more robust relationship between positive peaks and the muscele and posture signals, makes us conclude that in general the movements performed in our experiment are increasing subglottal pressures increasing the voice’s amplitude.

Show code producing the plot

sub <- tr_wd[ (tr_wd$trialindex>9) & (tr_wd$vocal_condition=='vocalize'),]
#make new dataset for faceting the plot
subpeak <- sub %>% 
  pivot_longer(cols = c("max_amp_time_around_move", "min_amp_time_around_move"),
               names_to = "peak_type",
               values_to = "time")

subpeak$peak_type <- recode(subpeak$peak_type,
                            max_amp_time_around_move = "positive peak",
                            min_amp_time_around_move = "negative peak")


densitydistance <- subpeak %>% 
  ggplot(aes(x = time, color = peak_type, linetype = movement_condition)) +
  geom_density(size = 0.75) +
  scale_linetype_manual(values = c("extension"="dotted", "external rotation"="dashed", "flexion"="dotdash", "internal rotation"="longdash", "no movement" = "solid")) +
  scale_colour_manual(values = c("red", "blue"), guide = "none") +
  xlab("Distance (in ms)") +
  ylab("Density") +
  labs(linetype="movement condition") +
  theme_cowplot() +
  theme(axis.text.x = element_text(angle = 45)) +
  background_grid() +
  facet_grid(~peak_type)

Figure S10: Density plot showing temporal distance (in milliseconds) of positive and negative peaks in amplitude relative to peak speed in wrist movement.