Data analysis

All analyses were performed in R (version 3.4.1; R Core Team, 2017).

During the data analysis, we first examined the general adaptation pattern that occurred over the course of the experiment in the syllables containing the central vowel /ɨ/. Next, we looked at individual spatial and temporal changes of vowel formants due to the applied perturbation. Finally, by investigating participants’ initial F1–F2 vowel space, we evaluated the potential influence of the surrounding sound categories /i/ and /u/ on the individual compensation strategies.

To examine average formant changes in participants’ production of the two syllables /dɨ/ and /gɨ/ across the four experimental phases, we fitted a generalized additive model (GAM; Hastie & Tibshirani, 1987). A GAM is a significant extension of a generalized linear regression model which allows the modelling of non-linear relationships between the depen-dent and independepen-dent variables (Wood, 2017a). Therefore, GAMs are much more flexible compared to linear regression models. The non-linear relationships are modelled via complex functions (smooths) which are constructed from ten basis functions (e.g., linear, quadratic, and cubic functions) with an adjustable number of basis dimensions. The number of basis dimensions is a number which indicates the upper limit of how com-plex the constructed function can be and is estimated directly from the data during the modelling process. That means that the usage of GAMs does not require from the researcher a predefined specification of a certain (non-linear) function as it is derived directly from the data. To prevent overfitting of the data, i.e., modelling of functions which are too com-plex and therefore might obscure any generalizable patterns in the data, GAMs are estimated using penalized likelihood estimation and cross-validation (cf. for details Wood, 2006). In the case of cross-cross-validation, several subsets of the complete data sample are created always excluding a single data point and the model is refitted to all of these subsets exam-ining how well it predicts the excluded data. One further advantage of GAMs is the possibility to include random effects into the model struc-ture to account for individual response variability across but also within speakers (cf. Baayen, Vasishth, Kliegl & Bates, 2017). To denote the inclusion of random effects in the fitted model, it is dubbed generalized

additive mixed model (GAMM). For a hands-on introduction to GAMMs with a focus on dynamic speech analysis see Sóskuthy (2017).

The GAMM offers three main advantages for analyzing the data from the current experiment. First, it is possible to analyze the data as a function of time which allows us to investigate the whole adaptation process rather than just its outcome. Secondly, the non-linearity of parameter smooths does not make any assumptions regarding the temporal or spatial charac-teristics of the adaptation process. Finally, the parameter smooths can be estimated including random effects which allows us to capture individual variability of the adaptation process.

Prior to building the GAMM model, participants’ raw formant fre-quencies were normalized by subtracting each participants’ mean formant frequency produced during the baseline phase for the respective syllable (/dɨ/ or /gɨ/). This was done to exclude participant-specific differences regarding their absolute formant magnitudes (e.g., due to gender differences). By means of this normalization, the average F1 and F2 values for /dɨ/ and /gɨ/ were set at zero for the baseline phase.

Subsequently, using the mgcv package (Wood, 2017b) we fitted one GAMM model for each formant (F1 and F2) with normalized frequencies averaged across all participants and all experimental trials as dependent variable. The data of the unperturbed syllables /di/ and /gu/, which were uttered by participants only during the baseline phase, were not included in the resulting GAMMs. All GAMM models were evaluated, interpreted, and visualized by means of the itsadug package by van Rij, Wieling, Baayen & van Rijn (2017).

In the model structure, we included random factor smooths with an intercept split for the perturbation direction (upward vs. downward) in order to assess (potentially non-linear) individual compensation mag-nitude differences over the course of the experiment. The model also included a fixed effect which assessed the ‘constant’ effect of the pertur-bation direction independently from the temporal variation. The resulting models explained 46.6 % and 66.9 % of the variance in the F1 and F2 data, respectively. In comparison, the model which did not include the random smooths (participant-specific temporal variation) but only random intercepts and random slopes explained only 31.2 % of the variance in the F2 data. Maybe somewhat surprisingly, the inclusion of the phase

number (shift 1, shift 2, and shift 3) as an interaction with the perturba-tion direcperturba-tion did not significantly improve the model fit. We also refitted the F2 model including an interaction between the perturbation direction (upward vs. downward) and the experimental group (A vs. B) which also did not improve the fit. In both cases, the goodness of fit was assessed by the Akaike Information Criterion (AIC; Akaike 1974).

Following the suggestion in Baayen, van Rij, de Cat & Wood (2016), the fitted models were investigated for the presence of autocorrelation in their residuals. Autocorrelation in the present study represents the correla-tion between the formant frequencies produced by one participant on two consecutive experimental trials. The higher the autocorrelation value is the less amount of information is contributed for the statistical model by each additional experimental trial. Ignoring this issue might result in overcon-fident estimates of the standard errors, confidence intervals, and p-values.

The amount of autocorrelation at lag 1 was relatively moderate in the present data with 0.2 for F1 and 0.17 for F2. The effect of autocorrela-tion was practically reduced to zero by incorporating AR(1) error models in the specification of the fitted GAMM models. The corrected models explained 23.1  % and 63.4  % of the variance in the F1 and F2 data, respectively. The dropped percentages of the explained variance are due to the refitted models taking into account the autocorrelation which makes their prediction about actual frequency values worse. This is especially true for the F1 model which is an indication that much of the variance in the initial model can be explained by autocorrelated errors rather than by the specified model parameters such as the direction of the applied pertur-bation. Visual model inspection revealed that the residuals of the adjusted GAMMs followed a normal distribution for F1 and F2 data.

To examine individual spatial and temporal differences of the adapta-tion process, we extracted F2 curves estimated for each participant by the GAMM model described in the above paragraphs.

In order to evaluate whether the occurrence of certain individual com-pensation patterns was induced by sound categories surrounding the perturbed vowel, we investigated participants’ F1–F2 space using their baseline phase production. For this purpose, we fitted two linear-mixed models using the lme4 package (Bates, Mächler, Bolker, and Walker, 2015).

One model was fitted for each of the two average formant frequencies (F1

and F2) that were produced by participants in the syllables /di/, /dɨ/, /gɨ/, and /gu/ during the baseline phase. The model structure included the pro-duced syllable and the interaction between the syllable and gender as fixed effects and the formant frequency as dependent variable. Furthermore, both models included an interaction between the syllable and the compen-satory pattern observed for each participant (cf. section 3.2 for a detailed discussion of individual compensation patterns). Random intercepts were modeled for each participant as well as random slopes for each produced syllable.

Visual model inspection revealed that the residuals of the chosen models followed a normal distribution for F1 and F2 data. P-values were obtained with the lmerTest package by Kuznetsova, Brockhoff, and Bojesen-Christensen (2016).

3. Results