1 Overview
This script contains the supplementary materials for the manuscript ‘Exploring the social meaning of the ’leader-lagger’ vowels in New Zealand English’. The pre-registration for this analysis can be accessed here.
The materials have the following structure:
Section 2 loads the R packages and data frames used in the analysis
Section 3 presents transcripts of the audio stimuli used in the pairwise similarity rating task, and a stock take of the monophthongs present in each stimulus.
Section 4 contains the code for the results and figures reported in the manuscript, specifically:
Fitting the reported MDS analysis in Section 4.2
Fitting a PCA to connecting MDS dimensions to participant labels in Section 4.3
Predicting MDS dimensions in Section 4.4
Section 6 presents the pre-registered correlations between the MDS dimensions and independent variables.
Finally, Section 7 discusses our data filtering decisions and presents the same analyses applied in Section 4 to the experiment data filtered according to the pre-registration.
If you would like to see how we scaled participant similarity ratings and created the similarity matrices used in this file, please refer to the file Matrix-generation-script.qmd in the github repository or view the file here.
2 Libraries and dataframes
The following code chunks load:
- Required libraries.
Click here to view code.
# Data manipulation
library(tidyverse)
# Set theme for visualisations
theme_set(theme_bw())
# Visualisations
library(ggcorrplot) # correlations
library(corrplot) # correlations
library(cowplot) # function plot_grid
library(ggpubr) # function ggarrange
library(rpart.plot) # visualise regression tree output
library(gratia) # visualise bam model
library(magick) # manipulate .png files
library(ggalt) # gg_encircle function
library(tidytext) # scale_y_reordered function
library(factoextra) # visualise label PCA output
# Statistical analyses
# MDS
library(smacof) # Apply MDS
library(rsample) # bootstrapping
library(vegan) # procrustes transformation of MDS
# Using development version of nzilbb.vowels,
# Any version after 0.3.1 will work.
# remotes::install_github('nzilbb/nzilbb_vowels')
library(nzilbb.vowels) # Assess MDS fit
# Decision trees
library(tidymodels) # fit decision trees
library(parsnip) # fit decision trees and random forests
library(randomForest) # run random forests
library(rpart) # run regression trees
library(vip) # assess importance of predictors in random forests
# Other
library(here) # localised file paths
library(knitr)
library(gt) # html tables when rendering
library(kableExtra) # html tables when rendering
library(grateful) # write package citations at end of documentTwo similarity matrices from the reported Free Classification task, converted to dissimilarity matrices for MDS analysis. One matrix uses the full, filtered, data. The second matrix uses a subset of the full data; when listeners used explicitly social labels. The values in both matrices represent the proportion of times each stimuli pair were placed in the same group of the times they could have been.
Original groups and labels for each group/speaker, for both the filtered and unfiltered full data set, and the social subset. The latter is used in the descriptive analysis of participant labels.
Click here to view code.
# Load filtered experimental results
label_df_full_filtered <- read_rds(
here(
"Data",
"df_full_anon_250228.rds"
)
)
paste0(
"Filtered participant count: ",
n_distinct(label_df_full_filtered$workerId)
)[1] "Filtered participant count: 116"Click here to view code.
# Load unfiltered experimental results
label_df_full_unfiltered <- read_rds(here(
"Data",
"df_full_unfiltered_anon_250228.rds"
))
paste0(
"Unfiltered participant count: ",
n_distinct(label_df_full_unfiltered$workerId)
)[1] "Unfiltered participant count: 127"Click here to view code.
# Load experimental results where listeners use social labels only
label_df_social <- read_rds(here(
"Data",
"df_social_anon_250228.rds"
))
# Bring personality and accent types into 'social factors'
# Bring 'NotSpeechNotSocial' and 'Unsure' (participants explictly being unsure) labels into single 'other' category
label_df_full_filtered <- label_df_full_filtered %>%
mutate(
label_category3 = case_when(
label_category2 == c("Personality") ~ "SocialFactors",
T ~ label_category3
),
label_category3 = case_when(
label_category3 == c("AccentTypes") ~ "SocialFactors",
T ~ label_category3
),
label_category3 = case_when(
label_category3 == "SocialFactors" ~ "Social",
label_category3 == "NotSpeechNotSocial" |
label_category3 == "Unsure" ~ "Other",
T ~ label_category3
)
)
label_df_full_unfiltered <- label_df_full_unfiltered %>%
mutate(
label_category3 = case_when(
label_category2 == c("Personality") ~ "SocialFactors",
T ~ label_category3
),
label_category3 = case_when(
label_category3 == c("AccentTypes") ~ "SocialFactors",
T ~ label_category3
),
label_category3 = case_when(
label_category3 == "SocialFactors" ~ "Social",
label_category3 == "NotSpeechNotSocial" |
label_category3 == "Unsure" ~ "Other",
T ~ label_category3
)
)
# Load social similarity matrix from experimental results.
# Matrix based on full data set (reported)
# Generated in Scripts/Matrix-generation-script.qmd
df_full <- read_rds(
here(
"Data",
"similarity_matrix_full_anon_250228.rds"
)
)
# Matrix based on subset of data (only when listeners used social categories)
# Also generated in Scripts/Matrix-generation-script.qmd
df_social <- read_rds(
here(
"Data",
"similarity_matrix_social_anon_250228.rds"
)
)
# Convert similarity matrices to dissimilarity matrices for MDS.
df_social <- sim2diss(df_social, method = 1)
df_full <- sim2diss(df_full, method = 1)- The similarity matrix from the Pairwise comparison task reported in Sheard et al. (In Prep), converted to a dissimilarity matrix for MDS analysis. Pairwise similarity ratings were scaled per participant, and the values in this matrix are the mean value, for each stimuli pair, of the scaled ratings. Please refer to this GitHub repository for data and scripts used to generated this matrix.
Click here to view code.
# Generated in separate GitHub repository
# Scaled pairwise ratings
matrix_scaled <- read_rds(
here(
"Data",
"PW_matrix_scaled_anon_250124.rds"
)
)
# Convert to dissimilarities for MDS.
df_PC <- sim2diss(matrix_scaled, method = "reverse")Additional data frames used in the analysis, specifically:
QB1 PCA scores
Measurements of articulation rate and pitch
Click here to view code.
# Load additional acoustic measures (PC loadings, GAMM intercepts, mean F1/F2, articulation rate, pitch)
QB1_scores_38 <- read_rds(here("Data",
"QB1_Scores_38_anon_250228.rds"))
# Calculate proportion of stimuli consisting of pause and duration
QB1_scores_38 <- QB1_scores_38 %>%
mutate(
pause_total = stimulus_duration - phonation_total,
phonation_prop = phonation_total / stimulus_duration * 100,
pause_prop = 100 - phonation_prop
) %>%
mutate(PC1_swapped = PC1 * (-1),
SpeakerID = as.character(SpeakerID))
# Mean normalised F1/F2 for whole monologue, outliers not removed
# For visualising example vowel spaces
QB1_vowels_mean <-
read_rds(here("Data", "QB1_vowels_mean.rds"))
# Load vowel counts for audio stimuli
VowelCount <- read_rds(here("Data",
"StimuliVowelFrequency_anon.rds"))
VowelCount <- VowelCount %>%
filter(!is.na(dress_E))2.1 Session information and additional functions
The code chunk below contains a custom function that creates a new column indicating whether a value in a specified column is above, below, or in between, a specified standard deviation, and applies a label to each of the three categories.
Click here to view code.
# Function to classify numeric variable as categorically low, middle or high in distribution
# based on SD (value)
sd_calculate <- function(orig_column, value, options) {
new_column <- case_when(
{{ orig_column }} > mean({{ orig_column }}, na.rm = TRUE) +
value * sd({{ orig_column }}, na.rm = TRUE) ~ options[3],
{{ orig_column }} < mean({{ orig_column }}, na.rm = TRUE) -
value * sd({{ orig_column }}, na.rm = TRUE) ~ options[1],
TRUE ~ options[2]
)
}The code chunk below contains a custom function that creates a correlogram in which the values in the bottom diagonal of the correlogram represent the correlation coefficients, and the values in the upper diagonal represent the p-values.
Click here to view code.
correlogram_function <- function(in_data) {
data_matrix <- as.matrix(in_data)
pmat <- cor.mtest(data_matrix, method = "spearman")
pvalues <- pmat$p
pvalues[upper.tri(pvalues)] <- NA
pvalues_long <- pvalues %>%
as.data.frame() %>%
tibble::rownames_to_column("Variable1") %>%
gather("Variable2", "value", -Variable1) %>%
mutate(
value = round(value, digit = 1),
value = case_when(Variable1 == Variable2 ~ NA, T ~ value)
)
test <- cor(data_matrix, method = "spearman")
corr_coeff <- test
corr_coeff[lower.tri(corr_coeff)] <- NA
corr_coeff_long <- corr_coeff %>%
as.data.frame() %>%
tibble::rownames_to_column("Variable1") %>%
gather("Variable2", "value", -Variable1) %>%
mutate(
value = round(value, digit = 1),
value = case_when(Variable1 == Variable2 ~ NA, T ~ value),
value = gsub("0.", ".", value)
)
correlogram_plot <- ggcorrplot(
test,
p.mat = pmat$p,
sig.level = 0.05,
insig = "pch",
pch = "",
method = "square",
outline = T,
type = "full",
# lab = TRUE,
# digits = 1,
show.diag = T,
ggtheme = ggplot2::theme_bw,
# lab_size = 2 #,
# tl.srt = 90
)
x_labs <- ggplot_build(
correlogram_plot
)$layout$panel_params[[1]]$x$get_labels()
correlogram_plot <- correlogram_plot +
scale_y_discrete(limits = rev) +
scale_x_discrete(position = "top", labels = x_labs) +
theme(axis.text.x = element_text(angle = 60, hjust = -0.05)) +
geom_text(
aes(
x = pvalues_long$Variable1,
y = pvalues_long$Variable2,
label = pvalues_long$value
),
size = 2,
fontface = "bold"
) +
geom_text(
aes(
x = corr_coeff_long$Variable1,
y = corr_coeff_long$Variable2,
label = corr_coeff_long$value
),
size = 2
)
correlogram_plot
}3 Full audio stimuli and vowel stock take
Table 1 displays the transcriptions for the 38 audio stimuli used in the online task. The table indicates whether each of the 10 NZE monophthongs considered in Brand et al. (2021) is present in the stimuli (1 = present, 0 = absent). The table also counts the total number of monophthongs represented in each stimulus (e.g., a 5 indicates 5 of the 10 monophthongs are present in a stimulus) and the total number of stimuli in which each monophthong is represented. We can see from the table that at least 6 of the 10 monophthongs are present in each stimuli, with a median of 8 per stimulus, and an average of 8.13 (Table 2 ). The vowels from the leader-lagger continuum are also more reliably represented than the vowels in the back-vowel configuration, largely due to the under-representation of the start vowel across the stimuli. We also note the strut F1 is loaded onto the leader-lagger continuum while strut F2 is loaded onto the back vowel configuration.
| speakerId | Stimulus | TRAP | DRESS | KIT | FLEECE | NURSE | STRUT | GOOSE | LOT | THOUGHT | START | Total |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | We tuned into the National Radio and the National Radio woman was fantastic she - she really deserves a - a bit of an award of how she handled it | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 7 |
| 2 | Everybody has gone through a lot of here - lot of problems here in Christchurch and we need to keep going | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 0 | 0 | 7 |
| 3 | Our older dog was a lot better actually because she's very deaf and quite blind so in fact that helped her a bit although obviously she knew something was up and she was outside, as well | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 8 |
| 4 | like deadliest catch I thought I was out in the middle of the sea with great huge waves and the noise was tremendous | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 7 |
| 5 | (the) only better thing I can think is the wonderful - relationships w- and people that I have come to know so much - so much better | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 6 |
| 6 | we sort of found room for everybody and everybody crashed the night um I don't think anybody slept we all slept in our clothes | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 0 | 8 |
| 7 | I was listening to National Radio and there was a wonderful um person there she was she was sort of making us all stay calm and things she was really good actually | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 9 |
| 8 | I couldn't get there cos there were too many cars and no-one everyone was just crawling like snails | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 10 |
| 9 | listening to the radio and all sitting round and talking so there were some good times and we had some laughs then as well | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 9 |
| 10 | they brought a bit of cutlery a plate because of course there was no water and we all just got together and um had this enormous meal and it was you know it was kind of nice um in a way | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 0 | 8 |
| 11 | I still have the people and the community um and one day those buildings'll be back in one form or another and in the meantime it's still a nice little community to be in | 1 | 0 | 1 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 6 |
| 12 | and i was trying really hard just to concentrate on the driving cos it was sort of the roads were starting to get quite sort of busy | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 10 |
| 13 | we've gained friends that we didn't know, who've come together under earthquake circumstances | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 8 |
| 14 | so at the end of it after two years we - we were in a very good position and we're both pretty happy with with what's happened for us | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 8 |
| 15 | and I went to go out in my car to go and find him and I got told it's not worth going out because there was so much traffic, there was so much damage | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 0 | 1 | 8 |
| 16 | I was with a group of friends and we had decided to leave at the end of the night and we walked out to shaking in the street | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 7 |
| 17 | there were so many people who came and helped from uh workmates to family to people who you don't even know | 1 | 1 | 0 | 1 | 1 | 0 | 1 | 1 | 0 | 0 | 6 |
| 18 | that was the th~ amazing thing about the earthquakes is that people actually asked if you were okay people wanted to know if things were okay . | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 0 | 1 | 7 |
| 19 | and in the end my dau~ daughter who is inc~ incredibly logical um just went to the phone book and tore out the back pages and we just slowly wor~ worked down the the list of things that we had to do . | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 9 |
| 20 | was quite hard to comprehend that this has happened as we're driving along in lovely weather and life is normal | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 9 |
| 21 | just immediately went into basically kiwi camping mode and I think that's what a lot of people did and it was just kind of second nature to lots of to to know how to do that | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 7 |
| 22 | for the September earthquake we were snuggled in bed and um we did all the things that we'd been taught as as children | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 7 |
| 23 | having said there were builders working on the building and they picked up all the children. The builders were great and they carried all the children into Latimer Square | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 9 |
| 24 | my initial thought was that is was something like an alien invasion because I had a sense that there was a- you know movement coming up from underneath us | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 0 | 8 |
| 25 | we found torches and a transistor radio, I turned on the transistor radio and there was um national program. on there was this really nice announcer and she just got information through from people | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 8 |
| 26 | I think I've got to a point now you just sorta look up check everything's alright nothing's damaged no one's hurt and just get . back on with it again | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 9 |
| 27 | I saw people walking past me so I was thinking I was thinking to myself well I'm gonna stop my car on the side of the road and walk because it's gonna be quicker . | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 8 |
| 28 | the afternoon was spent ahh gathering regrouping ourselves checking on . neighbours who were ahh shocked | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 0 | 1 | 8 |
| 29 | and here we have these four really nice council guys very big gentlemen . come over and here they are carrying this wedding dress across the road . struggling not to fall over . | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 9 |
| 30 | um . I was seeing my best friend that I hadn't seen for three years she was coming over from Honolulu so they wanted a real she wanted a real Sunday roast she'd been living in America | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 0 | 8 |
| 31 | I'm starting to feel . now some sense of . contentment and that it is my home again and it's got that feeling . of being home | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 7 |
| 32 | the next you know th~ few hours after that that horrible feeling of feeling um . sick and shakey and . um . that sort of feeling was awful but anyway life went on | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 9 |
| 33 | and all the cobble stones you know had to lift my feet because the cobblestones started . popping up . and I thought oh . that you know that might hurt if you got your foot . stuck or caught in that | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 9 |
| 34 | I always thought that even if I had had to walk all that way home . you could have trusted people cos people were so kind . and everyone said are you okay where are you walking to - | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 10 |
| 35 | so I managed to make nice coffee and some um . French toast and we found it c~ you know sorta felt quite guilty about having such delicious breakfast under such circumstances | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 10 |
| 36 | we've got a really lovely um . street where it's more of a village . a c~ and a community and people . really . stuck together and helped each other | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 0 | 8 |
| 37 | I stayed with my neighbours she cooked us food we all got together . it's funny when you don't know your neighbours that well . you bond pretty well - | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 0 | 8 |
| 38 | and and also just caring for people people . are a lot friendlier in Christchurch . we find people . are doing more for each other . . | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 10 |
| Total | 38 | 32 | 37 | 37 | 21 | 30 | 35 | 33 | 28 | 18 | NA |
Click here to view code.
VowelCount %>%
filter(!is.na(Total)) %>%
summarise(MeanCount = mean(Total),
MedianCount = median(Total)) %>%
kable() | MeanCount | MedianCount |
|---|---|
| 8.131579 | 8 |
4 Code for reported results
4.1 Participant labels
The code chunk below contains the code for Table 3, which displays the mean, median and ranges for both the number of groups made by participants, and the number of broad label categories used to describe these groups (Social, Speech or Other, see the third column of Table 1 in the manuscript). The subsequent chunk contains the code from Figure 2, which informs the discussion in Section 4.1 of the manuscript. Note that this analysis uses the unfiltered data.
Click here to view code.
groups <- label_df_full_unfiltered %>%
mutate(trial_index = as.factor(trial_index)) %>%
group_by(workerId, trial_index) %>%
# count number of groups made per participant and iteration
summarise(individual_label_count = n_distinct(ratings)) %>%
ungroup() %>%
summarise(
`Measure` = "Groups made",
Mean = mean(individual_label_count),
Median = median(individual_label_count),
Range = max(individual_label_count) - min(individual_label_count)
)
broad_cats <-label_df_full_unfiltered %>%
mutate(trial_index = as.factor(trial_index)) %>%
group_by(workerId, trial_index) %>%
summarise(
# count number of broad categories used per participant and iteration
broad_label_count = n_distinct(label_category3)
) %>% ungroup() %>%
summarise(
`Measure` = "Broad label categories used",
Mean = mean(broad_label_count),
Median = median(broad_label_count),
Range = max(broad_label_count) - min(broad_label_count)
)
# Create table with kable
groups %>%
rbind(broad_cats) %>%
kable() %>%
kable_styling(
bootstrap_options = "bordered",
position = "center",
font_size = 10
)| Measure | Mean | Median | Range |
|---|---|---|---|
| Groups made | 4.062992 | 4 | 12 |
| Broad label categories used | 1.671916 | 2 | 2 |
Click here to view code.
# Show number of different participants who used a certain label category
group_numbers <- label_df_full_unfiltered %>%
mutate(total_label_n = n()) %>%
group_by(label_category3) %>%
reframe(category_count = n(),
category_prop = category_count / total_label_n * 100) %>%
distinct() %>%
ungroup() %>%
arrange(category_prop) %>%
mutate(label_category3 = factor(label_category3, levels = label_category3)) %>%
ggplot(aes(x = label_category3, y = category_prop)) +
geom_segment(aes(xend = label_category3, yend = 0)) +
geom_point(size = 7, color = "orange") +
theme_bw(base_size = 12) +
labs(x = "Broad label category",
y = "Proportion of labels (%)",
title = "A") +
theme(
axis.text.x = element_text(
angle = 45,
vjust = 1,
hjust = 1
)) +
label_df_full_unfiltered %>%
mutate(trial_index = as.factor(trial_index)) %>%
group_by(workerId, trial_index) %>%
summarise(
individual_label_count = n_distinct(ratings),
broad_label_count = n_distinct(label_category3)
) %>%
ggplot(aes(y = individual_label_count, x = trial_index)) +
geom_boxplot() +
geom_jitter(alpha = 0.5) +
labs(x = "Iteration",
y = "Number of groups made",
title = "B") +
theme_bw(base_size = 12) +
label_df_full_unfiltered %>%
mutate(trial_index = as.factor(trial_index)) %>%
group_by(workerId, trial_index) %>%
summarise(
individual_label_count = n_distinct(ratings),
broad_label_count = n_distinct(label_category3)
) %>%
ggplot(aes(y = broad_label_count, x = trial_index)) +
geom_boxplot() +
geom_jitter(alpha = 0.5) +
labs(x = "Iteration",
y = "Number of broad label categories used",
title = "C") +
theme_bw(base_size = 12)
group_numbers
Click here to view code.
prop_label_fig <- label_df_full_unfiltered %>%
mutate(
label_category3 = case_when(
label_category3 == "SocialFactors" ~ "Social",
label_category3 == "NotSpeechNotSocial" |
label_category3 == "Unsure" ~ "Other",
T ~ label_category3
)
) %>%
group_by(label_category3, label_category2, label_category1) %>%
summarise(count = n_distinct(workerId)) %>%
group_by(label_category2) %>%
arrange(desc(count), .by_group = TRUE) %>%
mutate(label_category3 = factor(label_category3, levels = c("Speech", "Social", "Other"))) %>%
ggplot(aes(
x = count,
y = reorder_within(label_category1, count, label_category3)
)) +
geom_col(aes(fill = label_category2), colour = "black") +
facet_wrap( ~ label_category3, scales = "free_y", nrow = 1) +
scale_y_reordered() +
theme_bw(base_size = 10) +
theme(legend.position = "none") +
labs(y = "Label Category", x = "Number of participants")
# Save as png
ggsave(
prop_label_fig,
path = here("Figures", "Plots", "PNG"),
dpi = 300,
filename = "Figure4.png",
width = 3500,
height = 1500,
units = "px"
)
# Save as svg
ggsave(
prop_label_fig,
path = here("Figures", "Plots", "SVG"),
dpi = 600,
filename = "Figure4.svg",
units = "in",
height = 3.5
)
prop_label_fig
4.2 Fitting the reported MDS
Multi-dimensional scaling (MDS) is a dimension-reduction technique for similarity data used in multiple fields (e.g. psychology, marketing), with some implementation in linguistics (e.g., Clopper and Pisoni 2007; Clopper and Bradlow 2009). MDS reduces measures of (dis)similarity between pairs of objects to a smaller number of perceptual dimensions which, taken together, theoretically correspond to the cues or factors driving their perceived similarity (e.g., pitch, vowel patterns). We specifically implemented a monotone spline (M-spline) MDS, with five knots and third degree polynomials.1
MDS requires a dissimilarity matrix as input. For convenience, we first generated a square matrix which represents the similarity between each pair of speakers obtained from our experimental data. Each cell represents the number of times stimuli from two speakers were placed in the same group, in proportion to the number of times they could have been placed in the same group (stimuli were randomly distributed across three iterations of task, so not all participants grouped all the same stimuli). The similarity matrix was then converted to a dissimilarity matrix by subtracting each proportion from 1.
4.2.1 Determine number of dimensions
We now have a data frame for MDS. However, the number of dimensions we want must be specified in advance. A conventional approach is to use stress, a measure of fit (lower stress = better fit). However, because the addition of dimensions to an MDS analysis always reduces ‘stress’ (i.e., 3 dimensions will always have lower stress than 2 dimensions, and so on), directly choosing the number of dimensions from single stress values is a blunt tool and not considered best practice (see Mair, Borg, and Rusch 2016).
Here, we introduce a new method to inform the choice of the number of dimensions based on stress. This method is analogous to the permutation and bootstrapping approach developed by Viera (2012) and applied to vocalic co-variation by Wilson Black et al. (2023a). We implemented this method via a custom function in R in the nzilbb.vowels package. We want to ensure that we do not underestimate the dimensionality of the perceptual space. That is, we want to make sure we don’t set the dimensionality too low. In order to set a minimum number of dimensions we compare the stress reduction we would expect from permuted versions of our experimental data and for bootstrapped versions. We ensure that we have at least those dimensions which result in greater stress reduction than in permuted data.
Click here to view code.
set.seed(10)
full_test <- mds_test(
df_full,
n_boots = 100,
n_perms = 100,
test_dimensions = 5,
mds_type = "mspline",
spline_degree = 3,
spline_int_knots = 5
)
plot_mds_test(full_test)
Figure 3 suggests we need at least one dimension, with the black cross sitting somewhat higher than the null distribution for two dimensions. As previously noted, this test helps to convince us that we aren’t including too few dimensions.
4.2.2 Selecting best-fitting 2D analysis using random start values
The code chunk below uses our dissimilarity matrix to run 100 MDS analyses with random starts (M-spline MDS analyses, with five knots and third degree polynomials). The two-dimensional MDS analysis reported in the manuscript is the random start solution with the lowest stress value from these 100 runs. This step is motivated by Mair, Borg, and Rusch (2016), who note that there are many local minima for stress when fitting MDS.
Click here to view code.
set.seed(765)
fit_df_full <- NULL
for (i in 1:100)
fit_df_full[[i]] <- smacofSym(
df_full,
ndim = 2,
type = "mspline",
principal = T,
init = "random",
spline.degree = 3,
spline.intKnots = 5,
itmax = 2000
)
ind <- which.min(sapply(fit_df_full, function(x)
x$stress))
fit_df_full <- fit_df_full[[ind]]
fit_df_full
Call:
smacofSym(delta = df_full, ndim = 2, type = "mspline", init = "random",
principal = T, itmax = 2000, spline.degree = 3, spline.intKnots = 5)
Model: Symmetric SMACOF
Number of objects: 38
Stress-1 value: 0.305
Number of iterations: 410 4.2.3 Informal significance test assessing fit of final MDS
The informal significance test in Figure 4 gives an indication of the fit of the final MDS. The permtest function from the smacof package calculates the stress for 500 permuted iterations of the data frame, and the figure below visualises the distribution of the stress values across these iterations. The dotted line represents the lower 5% quantile, while the solid line represents the actual stress value of our final MDS analysis. The actual stress value is lower than the solid line, indicating that the stress of the actual MDS analysis is lower than over 95% of the analyses run on the permuted data.
Click here to view code.
set.seed(100)
nrep <- 500
res.perm_full <- permtest(
fit_df_full,
data = df_full,
method.dat = "maximum",
nrep = nrep,
verbose = FALSE
)
# permutation stress norm
mperm <- mean(res.perm_full$stressvec)
# lower 5% quantile (critical value)
perm5 <- quantile(res.perm_full$stressvec, probs = 0.05)
hist(
res.perm_full$stressvec,
xlim = c(0.10, 0.40),
xlab = "Stress Values",
main = "Permutation Histogram"
)
abline(v = perm5, lty = 2) # critical value (dashed)
abline(v = fit_df_full$stress)
The combination of Figure 4 and Figure 3 provides sufficient evidence for us to proceed with a two-dimensional MDS analysis.
4.2.4 Extracing speaker scores
We now extract the position of each speaker in the MDS space and join them with production information from the QuakeBox corpus.
Click here to view code.
# Extract the scores for each speaker.
conf_full <- fit_df_full$conf
dimensions_full <- as.data.frame(conf_full) %>%
rename(D1_full = V1, D2_full = V2) %>%
rownames_to_column(var = "Speaker")
# Join extracted scores with other data frames
QB1_scores_38 <- QB1_scores_38 %>%
right_join(dimensions_full, by = c("SpeakerID" = "Speaker"))
label_df_social2 <- label_df_social %>%
mutate(SpeakerID = as.character(SpeakerID)) %>%
right_join(dimensions_full, by = c("SpeakerID" = "Speaker"))The MDS scores for each speaker can be mapped in the two specified dimensions (see Figure 5). Stimuli which are close to one another in the 2D space are perceived to sound more similar to one another than those which are further apart.
Click here to view code.
QB1_scores_38 %>%
ggplot(aes(x = D1_full, y = D2_full)) +
geom_label(aes(label = SpeakerID)) +
theme_bw() +
labs(y = "Dimension 2 (D2)",
x = "Dimension 1 (D1)") +
coord_fixed() 
4.2.5 Rotating the MDS space to align with the airwise Rating space
The only thing which is important about this MDS space is the relative distance between the speakers. The specific dimensions of the space are arbitrary. While it is possible that stimuli D1 or D2 scores would directly align with a given feature in production, or with a different MDS space generated from a other data, this is not guaranteed. As we would like to compare the results of this analysis to those of Sheard et al. (In Prep), we can ensure that the MDS spaces are maximally comparable by rotating and scaling our MDS scores to maximally align the two spaces (i.e. we apply Procrustes rotation).
The code chunk below generates the same best fitting MDS analysis reported in Sheard et al. (In Prep)
Click here to view code.
set.seed(200)
fit_df_PC <- NULL
for (i in 1:100)
fit_df_PC[[i]] <- smacofSym(
df_PC,
ndim = 2,
type = "mspline",
principal = T,
init = "random",
spline.degree = 3,
spline.intKnots = 5,
itmax = 2000
)
ind <- which.min(sapply(fit_df_PC, function(x)
x$stress))
fit_df_PC <- fit_df_PC[[ind]]We extract the MDS scores and join with the QB data.
Click here to view code.
conf_PC_ID <- fit_df_PC$conf
dimensions_PW <- as.data.frame(conf_PC_ID) %>%
rename(D1_PW = V1, D2_PW = V2) %>%
rownames_to_column(var = "Speaker")
# Join extracted scores with other data
QB1_scores_38 <- QB1_scores_38 %>%
right_join(dimensions_PW, by = c("SpeakerID" = "Speaker"))We now use the procrustes function from the vegan package to rotate the space from the Free Classification task to be maximally similar to the space form the Pairwise ratings task.
Click here to view code.
fit_FC_rotated <- procrustes(fit_df_PC$conf,
fit_df_full$conf,
scores = "sites")We now extract the rotated scores and join them to the QB data.
Click here to view code.
rotated_FC <- fit_FC_rotated$Yrot %>%
as_tibble(.name_repair = "unique") %>%
rename(Rotated_D1_FC = `...1`,
Rotated_D2_FC = `...2`)
FC_dimensions <- dimensions_full %>%
select(-Speaker)
QB1_scores_38 <- QB1_scores_38 %>%
bind_cols(rotated_FC)The next chunk visualises the transformation of the Free Classification space to align with the Pairwise Comparison space, which is effectively a reversal of the horizontal axis (Figure 6).
Click here to view code.
rotation_1 <- QB1_scores_38 %>%
ggplot(
aes(
x = D1_full,
y = D2_full,
xend = Rotated_D1_FC,
yend = Rotated_D2_FC,
)
) +
geom_point(size=2) +
geom_segment(arrow = arrow(length = unit(2, "mm"))) +
labs(x="Dimension 1", y = "Dimension 2") +
coord_fixed() +
theme_bw()+
theme(legend.position="none")
rotation_1
We now put the two spaces side by side in Figure 7 (Figure 5 in the manuscript).
Click here to view code.
PC_space <- QB1_scores_38 %>%
ggplot(aes(x = D1_PW, y = D2_PW)) +
geom_label(aes(label = SpeakerID), size = 3) +
theme_bw() +
labs(y = "Dimension 2 (D2)",
x = "Dimension 1 (D1)",
title = "PW space (Sheard et al.)",
tag = "A") +
coord_fixed()
FC_orig <- QB1_scores_38 %>%
ggplot(aes(x = D1_full, y = D2_full)) +
geom_label(aes(label = SpeakerID)) +
theme_bw() +
labs(y = "Dimension 2 (D2)",
x = "Dimension 1 (D1)",
title = "Original FC space",
tag = "A") +
coord_fixed()
FC_procrustes <- QB1_scores_38 %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_label(aes(label = SpeakerID), size = 3) +
theme_bw() +
labs(y = "Dimension 2 (D2)",
x = "Dimension 1 (D1)",
title = "Free Classification") +
coord_fixed()
PC_space + FC_procrustes + plot_annotation(tag_levels = "A")
Click here to view code.
ggsave(
path = here("Figures", "Plots", "PNG"),
dpi = 300,
filename = "Figure5.png",
width = 2500,
height = 1500,
units = "px"
)
ggsave(
path = here("Figures", "Plots", "SVG"),
dpi = 300,
filename = "Figure5.svg",
# width = 5.5,
height = 3.5,
units = "in"
)4.2.6 Testing correlations between D1 and D2 scores from the Pairwise Rating and Free Classification spaces
The code chunk below tests a Pearson’s correlation between (i) D1 scores from the Pairwise Rating MDS space in Figure 7 A and rotated D1 scores from the Free Classification space in Figure 7 B and (ii) D2 scores from the Pairwise Rating MDS space in Figure 7 A and rotated D2 scores from the Free Classification space in Figure 7 B.
The correlation tests reveal a robust positive correlation between scores from both dimensions. From this, we can conclude that the Procrustes rotation has been effective in maximizing the similarity between the two spaces. We can also infer that the perceptual similarity spaces from the two tasks are indeed similar despite the different experimental mechanics. Section Section 5 explores the similarities between the two perceptual spaces in more detail.
Click here to view code.
cor.test(QB1_scores_38$D1_PW, QB1_scores_38$Rotated_D1_FC, method = "pearson")
Pearson's product-moment correlation
data: QB1_scores_38$D1_PW and QB1_scores_38$Rotated_D1_FC
t = 8.3246, df = 36, p-value = 6.542e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6636792 0.8980306
sample estimates:
cor
0.8112434 Click here to view code.
cor.test(QB1_scores_38$D2_PW, QB1_scores_38$Rotated_D2_FC, method = "pearson")
Pearson's product-moment correlation
data: QB1_scores_38$D2_PW and QB1_scores_38$Rotated_D2_FC
t = 6.4734, df = 36, p-value = 1.627e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5404432 0.8530833
sample estimates:
cor
0.7334141 4.3 PCA to determine relationships between MDS dimensions and labels
We now have our rotated perceptual similarity space, but what does it mean? We begin by exploring whether each dimension is interpretable. MDS as a technique is designed for similarity data, but it doesn’t tell you what the drivers of the similarity are. It’s left to researchers to identify what the underlying structure might be. Fortunately for us, we have our group labels! So the question for us is: as dimension scores change, which label categories correspondingly increase or decrease?
Speakers used so many different labels that differentiating between them quickly gets complex. To see which label categories also increase or decrease as D1/D2 scores increase or decrease, we implemented Principal Component Analysis (PCA). In this PCA the input data were D1/D2 scores and the proportion each social Level 1 label category (the most specific) was used to describe each stimuli (categories used by <10% of participants were excluded). To give an idea of what the data looks like, each row in Table 4 represents six of the 38 stimuli. Each column corresponds to a Level 1 social label categories, with each value representing the proportion that a given label category was used to describe that stimulus.
Click here to view code.
# I filtered out labels used by less than 10% of participants
# Then with the remaining labels calculated the proportion of labels for each participant that were from a particular category (e.g., 25% of their social labels referred to 'low SES')
label_test_full <- label_df_full_filtered %>%
mutate(SpeakerID = as.character(SpeakerID)) %>%
group_by(label_category1) %>%
mutate(nParticipant = n_distinct(workerId)) %>%
filter(nParticipant >= 12) %>%
group_by(SpeakerID) %>%
mutate(ParticipantTotal = n_distinct(workerId)) %>%
group_by(SpeakerID, label_category1) %>%
summarise(label_count = n(),
label_prop = label_count / first(ParticipantTotal) * 100) %>%
ungroup() %>%
select(SpeakerID, label_category1, label_prop) %>%
distinct() %>%
pivot_wider(id_cols = SpeakerID,
names_from = label_category1,
values_from = label_prop)
# The 0 is for NAs where a category wasn't used at all for a speaker (i.e., 0% of their labels refer to a particular category)
label_test_full[is.na(label_test_full)] <- 0
label_names <- label_test_full %>%
select(-SpeakerID) %>%
colnames()
label_PCA_full <- label_test_full %>%
right_join(QB1_scores_38, "SpeakerID") %>%
ungroup() %>%
select(all_of(label_names),
Rotated_D1_FC,
Rotated_D2_FC)
speaker <- label_test_full %>%
right_join(QB1_scores_38, "SpeakerID") %>%
ungroup() %>%
select(SpeakerID)| Breathy | British | ClearArticulation | Confident | DistVowels | ExpressiveTone | FastSpeech | FlatTone | Friendly | HighPitch | HighSES | Lisp | LowPitch | LowSES | MediumPitch | MiddleAge | MiddleSES | NZAverage | NZProper | NZRegional | NZRural | Nasal | Old | Quiet | SlowerSpeech | StrongNZ | UnclearSpeech | Young | Creaky | DistConsonsants | NegativeEmotions | Rotated_D1_FC | Rotated_D2_FC |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1.05 | 1.05 | 8.42 | 4.21 | 3.16 | 2.11 | 8.42 | 3.16 | 3.16 | 1.05 | 5.26 | 1.05 | 5.26 | 2.11 | 3.16 | 4.21 | 2.11 | 8.42 | 4.21 | 1.05 | 2.11 | 3.16 | 1.05 | 3.16 | 3.16 | 4.21 | 2.11 | 8.42 | 0.00 | 0.00 | 0.00 | -0.43 | 0.11 |
| 2.11 | 4.21 | 4.21 | 1.05 | 2.11 | 2.11 | 10.53 | 1.05 | 2.11 | 3.16 | 2.11 | 4.21 | 1.05 | 3.16 | 3.16 | 6.32 | 3.16 | 5.26 | 0.00 | 1.05 | 1.05 | 7.37 | 6.32 | 3.16 | 2.11 | 4.21 | 2.11 | 1.05 | 5.26 | 4.21 | 1.05 | -0.16 | 0.46 |
| 0.00 | 2.13 | 11.70 | 4.26 | 2.13 | 2.13 | 9.57 | 5.32 | 2.13 | 1.06 | 2.13 | 0.00 | 4.26 | 1.06 | 2.13 | 3.19 | 6.38 | 11.70 | 1.06 | 2.13 | 0.00 | 1.06 | 4.26 | 5.32 | 2.13 | 1.06 | 1.06 | 4.26 | 2.13 | 4.26 | 0.00 | -0.45 | 0.00 |
| 1.01 | 2.02 | 13.13 | 4.04 | 1.01 | 5.05 | 6.06 | 4.04 | 0.00 | 0.00 | 4.04 | 3.03 | 10.10 | 0.00 | 3.03 | 2.02 | 4.04 | 6.06 | 4.04 | 2.02 | 0.00 | 2.02 | 11.11 | 1.01 | 2.02 | 1.01 | 0.00 | 1.01 | 3.03 | 4.04 | 0.00 | -0.04 | -0.01 |
| 3.16 | 4.21 | 8.42 | 1.05 | 3.16 | 3.16 | 0.00 | 3.16 | 2.11 | 3.16 | 6.32 | 2.11 | 8.42 | 0.00 | 0.00 | 3.16 | 3.16 | 3.16 | 4.21 | 0.00 | 0.00 | 1.05 | 7.37 | 9.47 | 10.53 | 1.05 | 1.05 | 2.11 | 2.11 | 3.16 | 0.00 | -0.14 | -0.57 |
| 0.96 | 0.00 | 10.58 | 1.92 | 1.92 | 1.92 | 0.96 | 4.81 | 0.00 | 0.00 | 5.77 | 3.85 | 8.65 | 0.96 | 0.96 | 0.96 | 2.88 | 9.62 | 3.85 | 0.00 | 0.00 | 1.92 | 4.81 | 5.77 | 15.38 | 1.92 | 3.85 | 2.88 | 0.96 | 0.96 | 0.96 | -0.25 | -0.54 |
4.3.1 Determining the number of Principal Components
While Principal Component Analysis does not require researchers to specify the number of dimensions in advance, we do need to decide how many principal components should be considered in the interpretation of its results. As such, we again use a permutation and bootstrapping procedure that compares the variance explained by each PC in bootstrapped versus randomised data (see Wilson Black et al. 2023a for more information on this method and how to read the charts in this section). Figure 8 shows that the variance explained by the principal components from PC3 onwards for the bootstrapped data (red, with the dots indicating values obtained from the full data set) do not explain more variance than would be explained by chance (blue). As with the MDS, we decide to stick with two Principal Components.
Click here to view code.
PCA_full <- prcomp(label_PCA_full, scale = T)
PCA_test_full <- pca_test(
label_PCA_full,
n = 1000,
variance_confint = 0.95,
loadings_confint = 0.9
)
plot_variance_explained(PCA_test_full)
Click here to view code.
# flip the first PC (pca_test)
flipped_pca_test <-
pc_flip(PCA_test_full, pc_no = 1, flip_var = "NZRural")
# flip the second PC (pca_test)
flipped_pca_test <-
pc_flip(flipped_pca_test, pc_no = 2, flip_var = "SlowerSpeech")
# flip the first PC (prcomp)
flipped_pca <- pc_flip(PCA_full, pc_no = 1, flip_var = "NZRural")
# flip the second PC (prcomp)
flipped_pca <-
pc_flip(flipped_pca, pc_no = 2, flip_var = "SlowerSpeech")4.3.2 Labels and dimensions loaded onto PC1:
Figure 9 depicts the variables loaded onto PC1 from the label PCA (Figure 6A in the manuscript). A ‘+’ indicates that a variable is positively loaded on the PC (as PC1 score increases, the value of the variable increases). A ‘-’ indicates that label category is negatively loaded on the PC (as PC1 score increases, value of that variable decreases). The figure also represents the estimated randomised/permuted (blue) and bootstrapped (red) distribution of the loadings for each variable. If an index loading (+/- sign) falls above the 90% confidence band for the null distribution (i.e., above the blue line), then we consider it meaningfully captured by the PC.
Click here to view code.
PC1_loadings <- plot_loadings(flipped_pca_test, pc_no = 1) +
theme_bw(base_size = 12) +
labs(x = "Labels associated with D1") +
theme(
plot.title = element_blank(),
plot.subtitle = element_blank())
PC1_loadings <- PC1_loadings + theme(legend.position = 'none')
PC1_loadings
D1 (Rotated_D1_FC, with FC referring to the Free Classification task) is loaded onto PC1, along with the labels NZ Rural (accent), Strong NZ (accent), Low, Middle and High Socioeconomic status (SES), and Clear articulation, Quiet, and Creaky as the most strongly loaded speech-related labels. Young, Distinct Vowels and NZ Average also meet our 90% confidence band level threshold. The loadings indicate that, as stimuli D1 scores increase, use of High and Middle SES labels decreases and use of Low SES labels and rural labels increases. Dimension 1 from the MDS therefore appears to primarily capture perceived social class and accent broadness.
4.3.3 Labels and dimensions loaded onto PC2
Figure 10 is the same as Figure 9 except it depicts the variables loaded onto PC2 from the label PCA (Figure 6B in the manuscript).
Click here to view code.
PC2_loadings <- plot_loadings(flipped_pca_test, pc_no = 2) +
theme_bw(base_size = 12) +
labs(x = "Labels associated with D2") +
theme(plot.title = element_blank(),
plot.subtitle = element_blank())
PC2_loadings <- PC2_loadings + theme(legend.position = 'none')
PC2_loadings
D2 (Rotated_D2_FC) is loaded onto PC2, with the labels Fast Speech, Slower Speech, Low, Medium and High pitch, and Old, NZ Proper (accent) and Flat tone also meeting our 90% confidence band level threshold.2 As D2 scores increase, uses of Slow and Low Pitch labels increase. As D2 scores decrease, uses of Fast and High/Middle pitch labels increase. Dimension 2 from the MDS analysis therefore appears to capture perceived speed and pitch.
We extract scores from the label PCA and join with the QB data.
Click here to view code.
# Extract the scores and join with the QB information
coordinates_full <- get_pca_ind(flipped_pca)
coordinates_full <- as_tibble(coordinates_full$coord) %>%
select(Dim.1, Dim.2) %>%
cbind(speaker) %>%
rename(SES_Dim = Dim.1, Speed_Dim = Dim.2) %>%
left_join(QB1_scores_38, by = "SpeakerID")4.3.4 Connecting PCA to the MDS dimensions
What we found in the label PCA was a split between primarily social labels and labels primarily related to speed and pitch. This split also aligns with the two rotated dimensions from the MDS analysis, where Dimension 1 differentiates between the perceptually low SES and perceptually high SES speakers and Dimension 2 differentiates between the perceptually high/fast and perceptually slow/low speakers. Figure 11 shows how the label PCA translates to our two dimensions. Perceptually low SES speakers are encircled in yellow, perceptually high SES in purple. Perceptually fast/high speakers are encircled in red, perceptually low/slow speakers in blue.
Click here to view code.
coordinates_full %>%
mutate(
Speed_Dim = as.numeric(Speed_Dim),
SES_Dim = as.numeric(SES_Dim),
SES_Dim_Cat = sd_calculate(SES_Dim, 1, c("HighSES", "Middle", "LowSES")),
SES_Dim_Cat2 = sd_calculate(SES_Dim, 0.5, c("HighSES", "Middle", "LowSES")),
Speed_Dim_Cat = sd_calculate(Speed_Dim, 1, c("FastHigh", "Middle", "SlowLow")),
Speed_Dim_Cat2 = sd_calculate(Speed_Dim, 0.5, c("FastHigh", "Middle", "SlowLow"))
) %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_label(aes(label = SpeakerID), size = 3) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2),
alpha = 0.15
) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat != "Middle"),
aes(fill = Speed_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(fill = Speed_Dim_Cat2),
alpha = 0.15
) +
scale_fill_manual(values = c("#f0496a", "#a541f7", "#ffd117", "#09c9c3")) +
theme_bw() +
labs(
y = "Rotated Dimension 2 (D2)",
x = "Rotated Dimension 1 (D1)",
colour = "Perceptual Groups",
fill = "Perceptual Groups"
)
Figure 12 then combines the loadings from the label PCA with the rotated MDS space to isolate the specific speakers that are most associated with the perceptually high and low SES groups (Figure 6C in the manuscript) and the with the perceptually fast/high and slow/low groups (Figure 6D in the manuscript).
Click here to view code.
perceptual_broadness <- coordinates_full %>%
mutate(
Speed_Dim = as.numeric(Speed_Dim),
SES_Dim = as.numeric(SES_Dim),
SES_Dim_Cat = sd_calculate(SES_Dim, 1,
c("HighSES", "Middle", "LowSES")),
SES_Dim_Cat2 = sd_calculate(SES_Dim, 0.5,
c("HighSES", "Middle", "LowSES")),
Speed_Dim_Cat = sd_calculate(Speed_Dim, 1,
c("FastHigh", "Middle", "SlowLow")),
Speed_Dim_Cat2 = sd_calculate(Speed_Dim, 0.5,
c("FastHigh", "Middle", "SlowLow"))
) %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 == "Middle"),
shape = 4,
size = 3,
alpha = 0.25
) +
geom_point(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat, shape = SES_Dim_Cat),
size = 5,
alpha = 1
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2, shape = SES_Dim_Cat2),
size = 5,
alpha = 0.35
) +
annotate(
"text",
x = -0.4,
y = -0.15,
label = "paste(bold('-'))",
parse = TRUE,
colour = "black",
size = 16
) +
annotate(
"text",
x = 0.5,
y = 0,
label = "paste(bold('+'))",
parse = TRUE,
colour = "black",
size = 16
) +
scale_shape_manual(values = c(21, 22)) +
scale_fill_manual(values = c("#a541f7", "#ffd117")) +
scale_colour_manual(values = c("#a541f7", "#ffd117")) +
theme_bw() +
labs(
y = "Dimension 2 (D2)",
x = "Dimension 1 (D1)",
colour = "Perceptual Group",
fill = "Perceptual Group",
shape = "Perceptual Group"
) +
coord_fixed() +
theme(legend.position = "bottom")
perceptual_speed <- coordinates_full %>%
mutate(
Speed_Dim = as.numeric(Speed_Dim),
SES_Dim = as.numeric(SES_Dim),
SES_Dim_Cat = sd_calculate(SES_Dim, 1,
c("HighSES", "Middle", "LowSES")),
SES_Dim_Cat2 = sd_calculate(SES_Dim, 0.5,
c("HighSES", "Middle", "LowSES")),
Speed_Dim_Cat = sd_calculate(Speed_Dim, 1,
c("FastHigh", "Middle", "SlowLow")),
Speed_Dim_Cat2 = sd_calculate(Speed_Dim, 0.5,
c("FastHigh", "Middle", "SlowLow"))
) %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat != "Middle"),
aes(fill = Speed_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(fill = Speed_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat2 == "Middle"),
shape = 4,
size = 3,
alpha = 0.25
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat != "Middle"),
aes(fill = Speed_Dim_Cat, shape = Speed_Dim_Cat),
size = 5,
alpha = 1
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(fill = Speed_Dim_Cat2, shape = Speed_Dim_Cat2),
size = 5,
alpha = 0.35
) +
annotate(
"text",
x = -0.15,
y = -0.5,
label = "paste(bold('+'))",
parse = TRUE,
colour = "black",
size = 16
) +
annotate(
"text",
x = -0.15,
y = 0.4,
label = "paste(bold('-'))",
parse = TRUE,
colour = "black",
size = 16
) +
theme_bw() +
scale_shape_manual(values = c(23, 25)) +
scale_fill_manual(values = c("#f0496a", "#09c9c3")) +
scale_colour_manual(values = c("#f0496a", "#09c9c3")) +
labs(
y = "Dimension 2 (D2)",
x = "Dimension 1 (D1)",
colour = "Perceptual Group",
fill = "Perceptual Group",
shape = "Perceptual Group"
) +
coord_fixed() +
theme(legend.position = "bottom")
broadness_speed <- plot_grid(perceptual_broadness,perceptual_speed, nrow=1, labels = c("C","D"))
broadness_speed
Finally, the code chunk below combines Figure 9, Figure 10, and Figure 12 to make Figure 6 from the manuscript.
Click here to view code.
# Output perceptual groups as a plot.
ggsave(
plot = broadness_speed,
path = here("Figures",'Plots', "PNG"),
dpi = 400,
filename = "PerceptualGroups.png",
width = 2500,
height = 1250,
units="px"
)
combined <- plot_grid(shared_loadings, broadness_speed, ncol=1)
ggsave(
plot = combined,
path = here("Figures",'Plots', "PNG"),
dpi = 300,
filename = "Figure6.png",
width = 3000,
height = 2400,
unit="px"
)
ggsave(
plot = combined,
path = here("Figures",'Plots', "SVG"),
dpi = 400,
filename = "Figure6.svg",
width = 10,
height = 8
)4.4 Regression trees to investigate links between perceptual dimensions and acoustic factors
MDS analyses often use pairwise correlations to identify acoustic correlates of perceptual dimensions (e.g., Nolan, McDougall, and Hudson 2011; McDougall 2013). We fit regression trees rather than pairwise correlations as they allow us to look at how multiple variables predict our perceptual dimensions simultaneously. Regression trees partition data sets into smaller groups (‘nodes’) and then fit a model for each subgroup based on the specified predictors. The tree creates an if-else statement for the most important predictor at each node that partitions the data into further nodes, until specified thresholds are met. For each node in a tree, we can see the estimated dependent variable, the proportion and the number of observations represented in the node.
This approach was chosen over single linear regressions due to their capacity to capture potential non-linear relationships between independent and dependent variables, and interactions between independent variables. Furthermore, the if-else statements regression trees produce effectively provide ‘cutoffs’ that, in our case, are very useful for delimiting groups of speakers within the MDS space (i.e., speakers above a certain cutoff should be concentrated in a similar area of the MDS perceptual space).
We specifically fit two regression trees with stimuli rotated D1 and D2 scores as the dependent variables, with the following independent variables: speakers’ scores that measure their position in the leader-lagger continuum and back-vowel configuration (generated in the analysis in Hurring et al. Under review), and stimuli articulation rate and mean pitch measures. All independent variables were scaled to facilitate comparability.
To mitigate the high variance and instability of regression trees, we took a conservative approach to fitting each tree. To mitigate the tendency of regression trees to overfit data, the settings specified there had to be an increase of at least 0.05 to the cost/complexity parameter for a split to occur. The regression trees were fit using the parsnip package (Kuhn and Vaughn 2024) in R with the mode set to regression and the engine set to Rpart (Thernau, Atkinson, and Ripley 2023). We also implemented random forests to evaluate the importance of each predictor in predicting each dimension. The random forests were also fit using parsnip in R with the mode set to regression and the engine set to Ranger (Wright and Ziegler 2017). We fit 1000 trees with the same improvements to the cost/complexity parameter for a split to occur as the regression tree.
4.4.1 Predicting D1 with a regression tree
Figure 13 depicts the results of the regression tree predicting Dimension 1 fit in the code chunk below. Starting with all 38 stimuli, we can see that position in leader-lagger continuum is the first node in the tree. Speakers who are laggers in the leader-lagger continuum with scores below -0.26 (more negative scores = further behind in changes) are estimated to have lower D1 scores than speakers with scores above 0.26. As negative D1 scores are associated with higher SES labels, this points to laggers in change being perceived to sound higher SES than leaders.
Within the leaders, however, their scores are mediated by their articulation rate: it is the slower leaders who have the highest estimated D1 scores. In other words, the speakers who are perceived to sound the most broad and low SES are not only leaders in change, they are also slower. Leaders in change who are fast are perceived to sound both less broad than slow leaders. As such, the primary differentiation of perceived broadness is between laggers and slow leaders.
Click here to view code.
# We set up a data frame for regression analysis, scaling all relevant
# variables.
regression_df <- coordinates_full %>%
mutate(
LeaderLaggerScore = scale(PC1_swapped),
BackVowelScore = scale(PC2),
ArticRate = scale(articulation_rate),
MeanPitch2 = scale(pitch_cross_75_500),
SpeechRate = scale(speech_rate_update),
PhonProp = scale(phonation_prop),
PauseProp = scale(pause_prop),
PauseTotal = scale(pause_total),
CreakProp = scale(creak_prop),
Syllables = scale(syllable_total),
Speed_Dim = as.numeric(Speed_Dim),
SES_Dim = as.numeric(SES_Dim),
SES_Dim_Cat = sd_calculate(SES_Dim, 1, c("HighSES", "Middle", "LowSES")),
SES_Dim_Cat2 = sd_calculate(SES_Dim, 0.5, c("HighSES", "Middle", "LowSES")),
Speed_Dim_Cat = sd_calculate(Speed_Dim, 1, c("FastHigh", "Middle", "SlowLow")),
Speed_Dim_Cat2 = sd_calculate(Speed_Dim, 0.5, c("FastHigh", "Middle", "SlowLow"))
)
# Now fit the trees
set.seed(4)
tree_spec <-
# Set model type
decision_tree() %>%
# Set mode
set_mode("regression") %>%
# Specify engine
set_engine(
"rpart",
control = rpart.control(cp = 0.05),
method = "anova",
model = TRUE
)
d1_fit <- tree_spec %>%
fit(Rotated_D1_FC ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df)We output a plot of the D1 regression tree.
Click here to view code.
d1_fit %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#a541f7", "#F9A742", "#be4580", "#e66f70", "#ffd117"),
cex = 1,
xflip = F
)
The following code block saves the plot.
Click here to view code.
ppi <- 300
png(
here("Figures", 'Plots', "PNG", "D1RegressionTreeFull.png"),
width = 1900,
height = 1350,
res = ppi
)
d1_fit %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#a541f7", "#F9A742", "#be4580", "#e66f70", "#ffd117"),
cex = 1,
xflip = F
)
dev.off()Figure 14 then maps the cutoff values from the regression tree on the MDS space. We can see that the slow leaders are concentrated to the right of the space, with almost all within the perceptually low SES group. The laggers are concentrated to the left of the space, with most of the speakers within the perceptually high SES group also laggers. Fast laggers tend to be between the two groups, but concentrated more with the laggers than the slow leaders.
4.4.1.1 Evaluating the stability of D1 predictor importance with random forests
But how stable are these features as predictors of D1 scores? The next code chunk applies random forests to the same dependent and independent variables used in the regression tree. Random forest uses rand_forest() to define a model that creates a large number of decision trees, each independent of the others. The final prediction uses all predictions from the individual trees and combines them. We are then able to extract from the combined trees in both procedures the estimated ‘importance’ of each of our independent variables in predicting the dependent variable.
Click here to view code.
set.seed(2)
# Random forest settings
random_f <-
rand_forest(trees = 1000) %>%
set_engine("ranger",
importance = "permutation",
alpha = 0.05) %>%
set_mode("regression")
# Apply random forest to dependent and independent variables
d1_rf_fit <-
random_f %>%
fit(Rotated_D1_FC ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df)Figure 15 shows the ranked importance of the four variables in the random forest regression trees. The importance value is instead based on whether a variable has a positive effect on the predictor performance of the tree. A more positive score corresponds to increased importance in predicting the dependent variable (i.e., more frequent use to make key decisions within the decision trees and greater improvement to model fit when used). It should be noted, however that these measures are for capturing the importance of individual variables and are not designed to capture interactions (see Wright, Ziegler, and König 2016).
Click here to view code.
d1_rf_fit %>%
extract_fit_engine() %>%
vip() +
theme_bw(base_size = 12) +
labs(y = "Importance (Random Forest)")
The results of the random forests uphold the primary importance of articulation rate and the leader-lagger continuum in predicting perceptual broadness. Pitch and the back vowel configuration have the least importance of the four predictors.
To give an idea of how perceptual broadness relates to speakers’ positions in the leader-lagger continuum, Figure 16 illustrates the different vowel spaces of two slow speakers at opposite ends of the leader-lagger continuum and with contrasting D1 scores (vowel spaces are based on mean F1/F2 values from their entire QuakeBox monologues, not just from the stimuli in the experiment). Vowels not in the leader-lagger continuum are in black, and vowels in the leader-lagger continuum are in colour. Compared to the lagger in change (low D1 score = high perceived SES), the leader in change (high D1 score = low perceived SES) has high and front realisations of dress, trap, nurse and lot, and retracted realisations of kit and fleece. Both speakers are slow.
Click here to view code.
QB1_vowels_mean %>%
filter(Speaker %in% c("QB_NZ_F_529", "QB_NZ_F_870")) %>%
pivot_longer(cols = c(2:21), names_to = "Variable") %>%
separate(Variable, into = c("Formant", "Vowel")) %>%
pivot_wider(
id_cols = c("Speaker", "Vowel"),
names_from = "Formant",
values_from = "value"
) %>%
mutate(
LeaderLagger = case_when(
Speaker == "QB_NZ_F_870" ~ " Leader (perceptually low SES, speaker 31)",
T ~ "Lagger (perceptually high SES, speaker 28)"
)
) %>%
ggplot(aes(y = F1, x = F2, label = Vowel)) +
geom_label(data = . %>% filter(!Vowel %in% c(
"KIT", "TRAP", "DRESS", "NURSE", "FLEECE", "LOT"
)),
colour = "black") +
geom_label(data = . %>% filter(Vowel %in% c(
"KIT", "TRAP", "DRESS", "NURSE", "FLEECE", "LOT"
)), aes(colour = Vowel)) +
scale_y_reverse() +
scale_x_reverse(expand = expansion(0.2)) +
facet_wrap( ~ LeaderLagger) +
theme_bw(base_size = 10) +
theme(legend.position = "none")
4.4.1.2 Confirming relationship between D1 and Label PC1 scores
This section has so far shown that speakers’ rotated D1 scores systematically increase and decrease with the use of labels related to (perceived) high and low socioeconomic status (captured by Label PC1), and that speakers are differentiated along D1 by their speed and position in the leader-lagger continuum. Figure 17 confirms that the Label PC1 scores themselves are also predicted by the same variables, based on the same regression tree procedure in the previous section, with some minor adjustments to the cutoff values for the leader-lagger scores and speed. The change in Leader-Lagger score cutoff reflects three fast/high speakers (laggers in the previous tree, leaders here) who are outside the perceptually high SES speaker group, and a fourth lagger who was not in the more concentrated high-SES speaker group.
Click here to view code.
set.seed(4)
pc1_fit <- tree_spec %>%
fit(SES_Dim ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df)
pc1_fit %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#a541f7", "#F9A742", "#be4580", "#e66f70", "#ffd117"),
cex = 1,
xflip = F
)
Nonetheless, the overall interpretation of the space in Figure 18 remains stable. Figure 18 maps the cutoffs from the regression tree predicting D1 scores onto the perceptual space on the right, and the cutoffs from the regression tree predicting PC1 scores on the left. Across both spaces, the perceptually high SES speaker group consists predominantly of the same laggers, and the perceptually low SES group consists of the same leaders (mostly slow)
Click here to view code.
perceived_broadness +
theme(legend.position = "left") +
regression_df %>%
mutate(
PC1_tree = case_when(LeaderLaggerScore < -0.65 ~ "Lagger", T ~ "Leader"),
speed_tree = case_when(ArticRate < (-0.2) ~ "Slow", T ~ "Fast"),
PC1_speed_tree = case_when(
PC1_tree == "Lagger" ~ "Lagger",
T ~ interaction(PC1_tree, speed_tree)
),
) %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat, linetype = SES_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2, linetype = SES_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 == "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
alpha = 0.25
) +
scale_linetype_manual(values = c("solid", "dashed")) +
scale_fill_manual(values = c(c("#a541f7", "#ffd117"))) +
scale_colour_manual(values = c(c("#a541f7", "#F9A742", "#ffd117"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = NULL,
colour = "LL/Speed",
shape = "LL/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group",
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "none", # change legend title font size
legend.text = element_text(size = 6))
4.4.2 Predicting D2 with a regression tree
Figure 19 displays the regression tree predicting perceptual speed/pitch (D2). Here, the leader-lagger continuum was not selected by the regression tree and articulation rate is the first node in the tree. Slow speakers with rates below -0.11 (more negative scores = slower) are estimated to have lower D2 scores than speakers with scores above -0.11. As negative D2 scores are associated with slower speech labels, this points to slow speakers being perceived as such. Moreover, within the group of faster speakers, those with a higher mean pitch have a higher estimated D2 score than those with lower mean pitch. This points to fast and higher pitch speakers also being perceived as such, as high D2 scores are associated with higher pitch/faster labels.
Click here to view code.
set.seed(9)
d2_fit <- tree_spec %>%
fit(Rotated_D2_FC ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df)
d2_fit %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#09c9c3", "#8BDAB2", "#F9A742", "#F87233", "#f0496a"),
cex = 1,
xflip = T
)
The following code block saves the plot above.
Click here to view code.
ppi <- 300
png(
here("Figures",
"Plots",
"PNG",
"D2RegressionTreeFull.png"),
width = 1400,
height = 1400,
res = ppi
)
d2_fit %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#09c9c3", "#8BDAB2", "#F9A742", "#F87233", "#f0496a"),
cex = 1,
xflip = T
)
dev.off()Figure 20 then maps the cutoff values from the regression tree in Figure 19 on to the MDS space. We can see that the slow leaders are concentrated in the bottom of the space, with almost all within the perceptually slow/low group. The fast and high pitch speakers are concentrated in the top left of the space.
4.4.2.1 Evaluate stability of D2 predictor importance with random forests
The next code chunk applies the same random forests procedures from perceptual broadness to perceptual speed. Figure 21 shows the ranked importance of the same four variables in the random forest procedure. Articulation rate emerges as by far the most important predictor, with mean pitch having the next most important impact. The back-vowel configuration and leader-lagger scores do not emerge as importance predictors of perceptual speed and pitch.
Click here to view code.
set.seed(7)
d2_rf_fit <-
random_f %>%
fit(Rotated_D2_FC ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df)
d2_rf_fit %>%
extract_fit_engine() %>%
vip() +
theme_bw(base_size = 12) +
labs(y = "Importance (Random Forest)")
4.4.2.2 Confirming relationship between D2 and Label PC2 scores
This subsection has so far shown that speaker’s rotated D2 scores systematically increase and decrease with the use of labels related to (perceived) speed and pitch (captured by Label PC2), and that speakers are differentiated along D2 by their speed and pitch. The code chunks below confirm that the Label PC2 scores themselves are also predicted by the same variables, with some minor adjustments to the cutoff values. As with D1 and PC1, none of the speakers who comprise the perceptually high/fast or slow/low groups are affected and the interpretation of the space remains stable.
Click here to view code.
set.seed(9)
pc2_fit <- tree_spec %>%
fit(Speed_Dim ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df)
pc2_fit %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#09c9c3", "#8BDAB2", "#F9A742", "#F87233", "#f0496a"),
cex = 1,
xflip = T
)
When we plot these cutoffs in the MDS space in Figure 23 we can see that the perceptually slow/low group is still comprised of predominantly of slow speakers (with some low pitch speakers) and almost all fast and high speakers comprise the perceptually high/fast group (along with some fast but not high speakers).
4.5 Interpreting the full space by applying results of regression trees to MDS spaces
Figure 24 interprets the MDS space as a whole, drawing on the results of the regression trees predicting D1/D2 and the Label PCA (Figure 8C in the manuscript). We can discern four main overlaps in speaker production and listener perception:
Higher pitch and/or fast speakers are perceived as such, regardless of whether they are a leader or lagger, and are concentrated in the top right of the space in red.
Slower and lower pitch speakers are also perceived as such, and are concentrated in the bottom (left) of the space in blue.
Perceptually low SES speakers are predominantly slow leaders in change, and are concentrated to the right of the space in yellow.
Perceptually high SES speakers are predominantly laggers in change, and are concentrated to the left of the space in purple.
Click here to view code.
perceived_space <- regression_df %>%
mutate(
PC1_tree = case_when(LeaderLaggerScore < -0.26 ~ "Lagger", T ~ "Leader"),
speed_tree = case_when(ArticRate < (-0.11) ~ "Slow", T ~ "Fast"),
pitch_tree2 = case_when(MeanPitch2 < (0.15) ~ "Low", T ~ "High"),
main_groups4 = case_when(
speed_tree == "Fast" & Rotated_D2_FC > (0.25) ~ "High and/or fast",
PC1_tree == "Lagger" ~ "Laggers",
speed_tree == "Slow" &
Rotated_D2_FC < (-0.4) ~ "Slow and/or low",
speed_tree == "Slow" &
PC1_tree == "Leader" ~ "Leaders (Slow)",
speed_tree == "Fast" &
PC1_tree == "Leader" ~ "Leaders (Fast)",
T ~ interaction(pitch_tree2, speed_tree, PC1_tree)
)
) %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2, linetype = SES_Dim_Cat2),
alpha = 0.15,
size = 2
) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat != "Middle"),
aes(fill = Speed_Dim_Cat, linetype = Speed_Dim_Cat),
alpha = 0.35,
size = 2
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat, linetype = SES_Dim_Cat),
alpha = 0.35,
size = 2
) +
geom_point(aes(shape = main_groups4,
colour = main_groups4),
size = 6,
stroke = 2) +
labs(title = "Interpretation of perceptual space",
shape = "Main Groups",
fill = "Associated labels") +
#- fast, medium/high pitch,
#+ slow speech, low pitch, flat tone, old
scale_fill_manual(values = c(c(
"#f0496a", "#a541f7", "#ffd117", "#09c9c3"
))) +
scale_colour_manual(values = c(c(
"#f0496a", "#a541f7", "orange", "#ffd117", "#09c9c3"
))) +
scale_shape_manual(values = c(18, 16, 17, 15, 8, 21)) +
annotate(
"label",
label = c("Slow (+low) leaders + laggers"),
y = -0.55,
x = 0,
fill = "#3edec5"
) +
annotate(
"label",
label = c("'Slow','low','old'"),
y = -0.65,
x = 0,
fill = "#3edec5"
) +
annotate(
"label",
label = c("SLOW + LOW"),
y = -0.45,
x = 0,
colour = "#3edec5"
) +
annotate(
"label",
label = c("Slow or low"),
y = -0.1,
x = 0.5,
fill = "#F9A742"
) +
annotate(
"label",
label = c("'Rural','Low SES','Strong NZ'"),
y = -0.2,
x = 0.5,
fill = "#F9A742"
) +
annotate(
"label",
label = c("LEADERS"),
y = 0,
x = 0.5,
colour = "#F9A742"
) +
annotate(
"label",
label = c("Slow (+low)"),
y = -0.1,
x = -0.3,
fill = "#a541f7"
) +
annotate(
"label",
label = c("'High SES','Mid SES','Clear'"),
y = -0.2,
x = -0.3,
fill = "#a541f7"
) +
annotate(
"label",
label = c("LAGGERS"),
y = 0,
x = -0.3,
colour = "#a541f7"
) +
annotate(
"label",
label = c("High and/or fast leaders + laggers"),
y = 0.5,
x = -0.25,
fill = "#f0496a"
) +
annotate(
"label",
label = c("'High pitch', 'fast'"),
y = 0.4,
x = -0.25,
fill = "#f0496a"
) +
annotate(
"label",
label = c("HIGH + FAST"),
y = 0.6,
x = -0.25,
colour = "#f0496a"
) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
shape = "Production Groups",
colour = "Production Groups",
fill = "Perceptual Groups",
linetype = "Perceptual Groups"
) +
scale_linetype_manual(values = c("dashed", "solid", "dotted", "twodash")) +
coord_fixed() +
theme_bw() +
theme(
legend.position = "right",
#legend.key.size = unit(1, "cm"),
# change legend key size
#legend.key.height = unit(1, "cm"),
# change legend key height
#legend.key.width = unit(1, "cm"),
# change legend key width
legend.title = element_text(size = 12, face = "bold"),
# change legend title font size
legend.text = element_text(size = 12)
)
perceived_space
# save as png
ggsave(
perceived_space,
path = here("Figures", "Plots", "PNG"),
dpi = 300,
filename = "MainGroupsFullLabelsInterpretation.png",
width = 3200,
height = 1650,
units = "px"
)
The code chunk below combines Figure 24 with Figure 13, Figure 14, Figure 19 and Figure 20 for the manuscript (Figure 6).
Click here to view code.
# Load interpretation of full MDS space
maingroups <-
image_read(here("Figures", "Plots", "PNG",
'MainGroupsFullLabelsInterpretation.png'
)
)
# Load blank image (used to create balance in final image)
blank <-
image_read(here("Figures", "External", 'Blank.png'))
main_side <- image_crop(blank, geometry = "50x1200")
main_top <- image_crop(blank, geometry = "3400x200")
# Add blank space to top of interpretation of full space
to_append_main <- c(main_top, maingroups)
maingroups <- image_append(to_append_main, stack = T)
to_append_main2 <- c(main_side, maingroups)
maingroups <- image_append(to_append_main2, stack = F)
# Add a border
maingroups <-
image_border(maingroups, color = "#000000", geometry = "15x15")
# Load plot of output of regression tree predicting D1
D1RegressionTree <-
image_read(here("Figures", "Plots", "PNG",
'D1RegressionTreeFull.png'))
# Load plot of output of regression tree predicting D2
D2RegressionTree <-
image_read(here("Figures", "Plots", "PNG",
'D2RegressionTreeFull.png'))
# Rotate D2 tree
D2RegressionTree <- image_rotate(D2RegressionTree, degrees = 90)
#Load MDS space with cutoffs from D1 regression tree mapped
D1MDS <-
image_read(here("Figures", "Plots", "PNG",
'PerceivedBroadness.png'))
# Load MDS with cutoffs from D2 regression tree mapped
D2MDS <-
image_read(here("Figures", "Plots", "PNG",
'PerceivedSpeedNoLeg.png'))
D2MDS_leg <-
image_read(here("Figures", "Plots", "PNG",
'PerceivedSpeedLegend.png'))
# Combine D1 regression tree + mapped MDS space with blank space
D1_top <- image_crop(blank, geometry = "1900x350")
D1_bottom <- image_crop(blank, geometry = "1900x250")
to_append_D1 <- c(D1_top, D1RegressionTree, D1MDS, D1_bottom)
D1_appended <- image_append(to_append_D1, stack = T)
# Add borders
D1_appended <-
image_border(D1_appended, color = "#ffffff", geometry = "30x15")
D1_appended <-
image_border(D1_appended, color = "#000000", geometry = "15x15")
# Combine D2 regression tree + mapped MDS space with blank space
D2_top <- image_crop(blank, geometry = "3400x150")
D2_bottom <- image_crop(blank, geometry = "3400x50")
to_append_D2 <- c(D2MDS, D2RegressionTree)
D2_appended <- image_append(to_append_D2, stack = F)
to_append_D2_2 <- c(D2_appended, D2MDS_leg)
D2_appended_2 <- image_append(to_append_D2_2, stack = T)
to_append_D2_3 <- c(D2_top, D2_appended_2,D2_bottom)
D2_appended_3 <- image_append(to_append_D2_3, stack = T)
to_append_D2_4 <- c(main_side, D2_appended_3)
D2_appended_4 <- image_append(to_append_D2_4, stack = F)
D2_appended_4 <-
image_border(D2_appended_4, color = "#000000", geometry = "15x15")
# Connect D2 tree/MDS space image to interpretation of full space
to_append_MainD2 <- c(D2_appended_4, maingroups)
MainD2_appended <- image_append(to_append_MainD2, stack = T)
# Connect D2 tree/MDS space/full space image to D1 tree/MDS space
toappend_D1_MainD2 <- c(D1_appended, MainD2_appended)
D1_MainD2_appended <- image_append(toappend_D1_MainD2, stack = F)
# Add border
D1_MainD2_appended <-
image_border(D1_MainD2_appended,
color = "#000000",
geometry = "15x15")
# Annotate full image with additional text
D1_MainD2_appended <- D1_MainD2_appended %>%
image_annotate(
"A: Analysis of D1",
size = 150,
location = "+100+45",
color = "black",
degrees = 0,
boxcolor = "white"
) %>%
image_annotate(
"B: Analysis of D2",
size = 150,
location = "+2100+45",
color = "black",
degrees = 0,
boxcolor = "white"
) %>%
image_annotate(
"C: Combined interpretation of D1 and D2",
size = 150,
location = "+2100+1950",
color = "black",
degrees = 0,
boxcolor = "white"
)
D1_MainD2_appended2 <- D1_MainD2_appended %>%
image_annotate(
"Slow Leaders",
size = 70,
location = "+1450+1300",
color = "black",
degrees = 0,
boxcolor = "white"
) %>%
image_annotate(
"Fast Leaders",
size = 70,
location = "+750+1300",
color = "black",
degrees = 0,
boxcolor = "white"
) %>%
image_annotate(
"Laggers",
size = 70,
location = "+250+1300",
color = "black",
degrees = 0,
boxcolor = "white"
) %>%
image_annotate(
"Fast + high pitch",
size = 70,
location = "+4150+325",
color = "black",
degrees = 0,
boxcolor = "white"
) %>%
image_annotate(
"Fast + low pitch",
size = 70,
location = "+4150+900",
color = "black",
degrees = 0,
boxcolor = "white"
) %>%
image_annotate(
"Slow speakers",
size = 70,
location = "+4250+1400",
color = "black",
degrees = 0,
boxcolor = "white"
)
D1_MainD2_appended_resize <- image_scale(D1_MainD2_appended2, geometry="x2800")
image_write(
D1_MainD2_appended_resize,
path = here(
"Figures",
'Plots',
"PNG",
'Figure8.png'
),
format = "png",
density= "400",
quality=100
)
image_write(
D1_MainD2_appended2,
path = here(
"Figures",
'Plots',
"SVG",
'Figure8.pdf'
),
format = "pdf",
density= "400",
quality=100
)
png <- image_read(here(
"Figures",
'Plots',
"PNG",
'Figure8.png'
))
image_info(png)# A tibble: 1 × 7
format width height colorspace matte filesize density
<chr> <int> <int> <chr> <lgl> <int> <chr>
1 PNG 3916 2800 sRGB TRUE 1364303 157x157Click here to view code.
D1_MainD2_appended_test <- image_convert(png, format="svg")
image_write(
png,
path = here(
"Figures",
'Plots',
"SVG",
'Figure8.svg'
),
format = "svg",
density="400",
flatten=F,
quality=100
)
D1_MainD2_appended2
5 Comparing spaces across tasks
Some differences emerged in the results predicting Dimension 1 from this MDS space and the space reported in Sheard et al. (In Prep) based on pairwise comparison data from the same stimuli. This section contains code for Figures 9 and 10 from the manuscript, which compare perceptually high/low SES and perceptually fast/slow speaker groups (i.e., speakers whose Label PC scores are 0.5 SDs above/below the mean) in the reported Free Classification MDS space to those in the Pairwise Rating MDS space.
The code chunk below creates a dataframe for the figures, using cutoffs from the regression trees predicting rotated D1 and D2 scores from the Free Classification MDS space.
Click here to view code.
visualisation_df <- regression_df %>%
mutate(
PC1_tree = case_when(LeaderLaggerScore < -0.26 ~ 'Lagger', T ~ "Leader"),
speed_tree = case_when(ArticRate < (-0.2) ~ "Slow", T ~ "Fast"),
speed_tree2 = case_when(ArticRate < (-0.11) ~ "Slow", T ~ "Fast"),
speed_tree3 = case_when(ArticRate < (-0.27) ~ "Slow", T ~ "Fast"),
speed_tree_comb = case_when(
ArticRate < (-0.27) ~ "Slow",
ArticRate >= (-0.11) ~ "Fast",
T ~ "Middle"
),
PC1_tree_comb = case_when(
LeaderLaggerScore < (-0.26) ~ "Lagger",
LeaderLaggerScore >= (-0.065) ~ "Leader",
T ~ "Middle"
),
pitch_tree_comb = case_when(
MeanPitch2 < (0.15) ~ "Low",
MeanPitch2 >= (0.96) ~ "High",
T ~ "Middle"
),
pitch_tree2 = case_when(MeanPitch2 < (0.15) ~ "Low", T ~ "High"),
BVC_tree = case_when(BackVowelScore > 0 ~ "Positive", T ~ "Negative"),
pitch_speed_tree = case_when(
speed_tree2 == "Slow" ~ "Slow",
speed_tree2 != "Slow" &
pitch_tree2 == "High" ~ "FastHigh",
T ~ "FastLow"
),
pitch_speed_tree2 = interaction(pitch_tree2, speed_tree2),
pitch_speed_tree3 = case_when(
pitch_tree2 == "High" ~ "High",
pitch_tree2 != "High" &
speed_tree2 == "Fast" ~ "FastLow",
T ~ "SlowLow"
),
PC1_speed_tree = case_when(
PC1_tree == "Lagger" ~ "Lagger",
T ~ interaction(PC1_tree, speed_tree)
),
main_groups3 = case_when(
PC1_tree == "Lagger" &
SES_Dim_Cat2 == "HighSES" ~ "DistinctLaggers",
PC1_tree == "Leader" &
SES_Dim_Cat2 == "LowSES" ~ "DistinctLeaders",
T ~ "NonDistinct"
)
) %>%
separate(
Speaker,
into = c("Corpus", "NZ", "Gender", "Number"),
sep = "_",
remove = F
) Figure 25 compares the D1 scores from both spaces, using different point types/colours to identify speakers with different leader-lagger vowels (Figure 9 in the manuscript). The perceptually low SES leaders and perceptually high SES laggers are in yellow circles and purple squares, respectively (i.e., perceptually low SES laggers and high SES leaders are not indicated). We can see that the same perceptually distinct leaders are grouped together in both spaces. These leaders are also grouped apart from the perceptually distinct laggers, who are grouped together in both spaces.
Click here to view code.
visualisation_df %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat, linetype = SES_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2, linetype = SES_Dim_Cat2),
alpha = 0.35
) +
geom_point(
data = . %>% filter(main_groups3 != "NonDistinct"),
aes(shape = main_groups3,
colour = main_groups3),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(main_groups3 == "NonDistinct"),
size = 6,
stroke = 1,
shape = 4,
colour = "grey",
alpha = 0.25
) +
geom_text(data = . %>% filter(main_groups3 != "NonDistinct"),
aes(label = SpeakerID)) +
geom_text(
data = . %>% filter(main_groups3 == "NonDistinct"),
aes(label = SpeakerID),
alpha = 0.25
) +
scale_shape_manual(values = c(16, 15)) +
scale_fill_manual(values = c(c("#a541f7", "#ffd117"))) +
scale_colour_manual(values = c(c(
"#a541f7", "#ffd117"
))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "Pitch/Speed",
shape = "Pitch/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group",
title = "Free Classification Task",
tag = "A"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "bottom",
legend.text = element_text(size = 6)) +
guides(shape = "none",
colour = "none") +
visualisation_df %>%
ggplot(aes(x = D1_PW, y = D2_PW)) +
geom_point(
data = . %>% filter(main_groups3 != "NonDistinct"),
aes(shape = main_groups3,
colour = main_groups3),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(main_groups3 == "NonDistinct"),
size = 6,
stroke = 1,
shape = 4,
colour = "grey",
alpha = 0.25
) +
geom_text(data = . %>% filter(main_groups3 != "NonDistinct"),
aes(label = SpeakerID)) +
geom_text(
data = . %>% filter(main_groups3 == "NonDistinct"),
aes(label = SpeakerID),
alpha = 0.25
) +
scale_shape_manual(values = c(16, 15)) +
scale_colour_manual(values = c(c(
"#a541f7", "#ffd117"
))) +
labs(
x = "Dimension 1 (D1)",
y = "Dimension 2 (D2)",
colour = "Main Groups",
shape = "Main Groups",
fill = "Perceptual Group",
linetype = "Perceptual Group",
title = "Pairwise Comparison Task",
tag = "B"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "bottom",
legend.text = element_text(size = 6)) +
guides(fill = "none")
# Save figure
ggsave(
path = here("Figures","Plots","PNG"),
dpi = 300,
filename = "Figure9.png",
width = 2500,
height = 1200,
units = "px"
)
ggsave(
path = here("Figures","Plots","SVG"),
dpi = 300,
filename = "Figure9.svg",
width = 8,
height = 3.5,
units = "in"
)
Similarly, Figure 26 compares D2 from the Free Classification on the left and the Pairwise Rating task on the right. The figure uses different point types/colours to identify speakers with different speed and pitch (Figure 10 from the manuscript). Slow, but high-pitched, speakers are positioned (pink circles) differently across the spaces. In the free classification space they are grouped more with the other slow speakers (in blue) than with the other higher pitched speakers (dark pink squares). In the pairwise comparison space, by contrast, the high-pitched speakers are consistently grouped together (i.e. the pink circles are with the pink squares - except for speaker 2). While the fast and high pitch speakers, and slow and lower pitch speakers, are perceptually distinct in both tasks, the slow and high pitch speakers are less stable.
Other than the slow and higher pitch speakers, however, the same speaker groups emerge across both tasks. Slow speakers are grouped together in the bottom of both tasks, and fast and high speakers are grouped together in the top left of both tasks.
Click here to view code.
visualisation_df %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat != "Middle"),
aes(fill = Speed_Dim_Cat, linetype = Speed_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(fill = Speed_Dim_Cat2, linetype = Speed_Dim_Cat2),
alpha = 0.35
) +
geom_point(
data = . %>% filter(interaction(speed_tree2, pitch_tree2) != "Fast.Low"),
aes(
shape = interaction(speed_tree2, pitch_tree2),
colour = interaction(speed_tree2, pitch_tree2)
),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(interaction(speed_tree2, pitch_tree2) == "Fast.Low"),
aes(
shape = interaction(speed_tree2, pitch_tree2),
colour = interaction(speed_tree2, pitch_tree2)
),
size = 6,
stroke = 1,
alpha = 0.25
) +
geom_text(data = . %>% filter(interaction(speed_tree2, pitch_tree2) !=
"Fast.Low"),
aes(label = SpeakerID)) +
geom_text(
data = . %>% filter(interaction(speed_tree2, pitch_tree2) == "Fast.Low"),
aes(label = SpeakerID),
alpha = 0.25
) +
scale_shape_manual(values = c(18, 15, 8, 4)) +
scale_colour_manual(values = c(c("#f0496a", "pink", "#09c9c3", "grey"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "Pitch/Speed",
shape = "Pitch/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group",
title = "Free Classification Task",
tag = "A"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "bottom",
legend.text = element_text(size = 6)) +
guides(shape = "none",
colour = "none") +
visualisation_df %>%
ggplot(aes(x = D1_PW, y = D2_PW)) +
geom_point(
data = . %>% filter(interaction(speed_tree2, pitch_tree2) != "Fast.Low"),
aes(
shape = interaction(speed_tree2, pitch_tree2),
colour = interaction(speed_tree2, pitch_tree2)
),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(interaction(speed_tree2, pitch_tree2) == "Fast.Low"),
aes(
shape = interaction(speed_tree2, pitch_tree2),
colour = interaction(speed_tree2, pitch_tree2)
),
size = 6,
stroke = 1,
alpha = 0.25
) +
geom_text(data = . %>% filter(interaction(speed_tree2, pitch_tree2) !=
"Fast.Low"),
aes(label = SpeakerID)) +
geom_text(
data = . %>% filter(interaction(speed_tree2, pitch_tree2) == "Fast.Low"),
aes(label = SpeakerID),
alpha = 0.25
) +
scale_shape_manual(values = c(18, 15, 8, 4)) +
scale_colour_manual(values = c(c("#f0496a", "pink", "#09c9c3", "grey"))) +
labs(
x = "Dimension 1 (D1)",
y = "Dimension 2 (D2)",
colour = "Main Groups",
shape = "Main Groups",
fill = "Perceptual Group",
linetype = "Perceptual Group",
title = "Pairwise Comparison Task",
tag = "B"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "bottom",
legend.text = element_text(size = 6)) +
guides(fill = "none")
ggsave(
path = here("Figures", "Plots", "PNG"),
dpi = 300,
filename = "Figure10.png",
width = 2500,
height = 1200,
units = "px"
)
ggsave(
path = here("Figures", "Plots", "SVG"),
dpi = 300,
filename = "Figure10.svg",
width = 8,
height = 3.5,
units = "in"
)
6 Testing pre-registered correlations
This section presents the correlations that were pre-registered for this analysis. Figure 27 shows a summary of the correlations between Dimensions 1 and 2 and speakers’ leader-lagger scores, their BVC scores, their articulation rate, creak proportion, and their mean pitch.3
The results of the regression trees and random forests are upheld in that Leader Lagger scores significantly correlate with Dimension 1 scores, and Articulation rate significantly correlates with Dimension 2 scores. However, unlike the regression trees and random forest, these pairwise correlations do not give us a sense of how the different variables relate to one another in differentiating between speakers in the perceptual space (e.g., once fast speakers are accounted for, Leader-Lagger scores correlate more strongly with Dimension 1).
Click here to view code.
columns <-
c(
"LeaderLaggerScore",
"BackVowelScore",
"ArticRate",
"MeanPitch2",
"CreakProp",
"Rotated_D1_FC",
"Rotated_D2_FC"
)
correlogram1 <- regression_df %>%
as_tibble() %>%
select(all_of(columns))
correlogram_function(correlogram1)
7 The role of the leader-lagger continuum when listeners are explicitly using social labels
As MDS is a dimension-reduction technique, it is currently unclear whether the differentiation between laggers and slow leaders in the perceptual space is a byproduct of the space condensing the different approaches to grouping and labelling speakers (i.e., speech versus social), or a meaningful perceptual distinction.
To investigate this, we apply the analyses as detailed in Section 4.2 and Section 4.3 to the subset of the data in which listeners specifically use social labels. This subsection:
Applies MDS to the social data subset
- Applies the same dimension selection and random-start procedure in the MDS
Rotates the best-fitting MDS to align with the space reported for the full data set
Applies a Label PCA to the social data subset
Predicts the rotated D1 and D2 scores for the social data subset with regression trees, with the same independent variables as
- It also explores D2 further in response to inconsistent results across the regression trees
7.1 Determine number of dimensions for the social data subset
Like Figure 3, Figure 28 suggests we need at least one dimension, with the black cross sitting somewhat higher than the null distribution for two dimensions. As Figure 28 indicates that two dimensions is not too few for this dataset, we proceed with two dimensions for the analysis in this section.
Click here to view code.
set.seed(10)
social_test <- mds_test(
df_social,
n_boots = 100,
n_perms = 100,
test_dimensions = 5,
mds_type = "mspline",
spline_degree = 3,
spline_int_knots = 5
)
plot_mds_test(social_test) +
labs(y = "Reduction in Stress",
x = "Number of Dimensions in MDS Analysis",
colour = "Test Distribution") +
scale_color_manual(values = c("#a788ef", "#ec8e9e"))
7.5 Regression trees to investigate links between perceptual dimensions and acoustic factors (social subset)
7.5.1 Predicting social D1 with a regression tree
The following code blocks does some prepatory data wrangling.
Click here to view code.
# variable scaling.
regression_df_social <- coordinates_social %>%
mutate(
LeaderLaggerScore = scale(PC1_swapped),
BackVowelScore = scale(PC2),
ArticRate = scale(articulation_rate),
MeanPitch2 = scale(pitch_cross_75_500),
SpeechRate = scale(speech_rate_update),
PhonProp = scale(phonation_prop),
PauseProp = scale(pause_prop),
PauseTotal = scale(pause_total),
PhonTotal = scale(phonation_total),
CreakProp = scale(creak_prop),
Syllables = scale(syllable_total),
Age_Dim = as.numeric(Age_Dim),
SES_Dim = as.numeric(SES_Dim),
SES_Dim_Cat = sd_calculate(SES_Dim, 1, c("LowSES", "Middle", "HighSES")),
SES_Dim_Cat2 = sd_calculate(SES_Dim, 0.5, c("LowSES", "Middle", "HighSES")),
Age_Dim_Cat = sd_calculate(Age_Dim, 1, c("Old", "Middle", "YoungFriendly")),
Age_Dim_Cat2 = sd_calculate(Age_Dim, 0.5, c("Old", "Middle", "YoungFriendly")),
)
# Set some cutoffs for visualisation.
pitch_cutoff <- (-0.15)
artic_cutoff <- (-0.5)
artic_cutoff2 <- (-0.11)
pause_prop_cutoff <- (0.31)
pause_time_cutoff <- (-0.67)
PC1_cutoff <- (-0.036)
PC1_cutoff2 <- (-0.26)
PC2_cutoff <- (0.74)
PC2_cutoff2 <- (-0.63)
# set up a data frame for visualisation purposes.
regression_df_social_vis <- regression_df_social %>%
mutate(
Age_Dim = as.numeric(Age_Dim),
SES_Dim = as.numeric(SES_Dim),
SES_Dim_Cat = sd_calculate(SES_Dim, 1, c("LowSES", "Middle", "HighSES")),
SES_Dim_Cat2 = sd_calculate(SES_Dim, 0.5, c("LowSES", "Middle", "HighSES")),
Age_Dim_Cat = sd_calculate(Age_Dim, 1, c("Old", "Middle", "YoungFriendly")),
Age_Dim_Cat2 = sd_calculate(Age_Dim, 0.5, c("Old", "Middle", "YoungFriendly")),
artic_tree = case_when(ArticRate < artic_cutoff ~ "Slow", T ~ "Fast"),
artic_tree2 = case_when(
ArticRate < artic_cutoff ~ "Slow",
ArticRate > artic_cutoff2 ~ "Fast",
T ~ "Middle"
),
pitch_tree = case_when(MeanPitch2 < pitch_cutoff ~ "Low", T ~ "High"),
pause_prop_tree = case_when(
PauseProp < pause_prop_cutoff ~ "LowPauseProp",
T ~ "HighPauseProp"
),
BVC_tree = case_when(BackVowelScore < PC2_cutoff ~ "LowBVC", T ~ "HighBVC"),
PC1_tree = case_when(LeaderLaggerScore < PC1_cutoff ~ "Lagger", T ~ "Leader"),
PC1_tree2 = case_when(
LeaderLaggerScore < PC1_cutoff2 ~ "Lagger",
LeaderLaggerScore > PC1_cutoff ~ "Leader",
T ~ "Middle"
),
BVC_tree2 = case_when(BackVowelScore < PC2_cutoff2 ~ "LowBVC2", T ~ "HighBVC2"),
BVC_tree3 = case_when(
BackVowelScore > PC2_cutoff ~ "HighBVC2",
BackVowelScore < PC2_cutoff2 ~ "LowBVC2",
T ~ "Middle"
),
speed_pause = case_when(
artic_tree == "Slow" ~ "Slow",
pause_prop_tree == "HighPauseProp" ~ "HighPauseProp",
T ~ artic_tree
),
PC1_speed_pause = case_when(
pause_prop_tree == "HighPauseProp" ~ "HighPauseProp",
artic_tree == "Slow" ~ "Slow",
T ~ interaction(pause_prop_tree, PC1_tree)
),
pitch_speed = interaction(artic_tree, pitch_tree),
pause_prop_time_tree = case_when(artic_tree == "Slow" ~ "Slow", T ~
interaction(pause_prop_tree))
)We now fit the regression tree.
Click here to view code.
set.seed(9)
fit_D1_social <-
tree_spec %>% fit(
Rotated_D1_FC_social ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df_social
)Figure 31 depicts the results of the regression tree predicting social D1 fit in the code chunk below. Starting with all 38 stimuli, we can see that position in leader-lagger continuum is, again, the first node in the tree. Speakers who are laggers in the leader-lagger continuum with scores below -0.32 (more negative scores = further behind in changes) are estimated to have lower D1 scores than speakers with scores above -0.32.
Within the leaders, their scores are, again, mediated by speed rate: the slower leaders have the highest estimated D1 scores. Slow leaders are the speakers who are perceived to sound the most broad and low SES. We therefore have additional evidence that the effect of the leader-lagger continuum on perceived socioeconomic status is mediated by the broader speech context, even when listeners are explicitly using social categories to group speakers.
Click here to view code.
fit_D1_social %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#a541f7", "#F9A742", "#be4580", "#e66f70", "#ffd117"),
cex = 1
)
When we compare the full and social MDS spaces in Figure 32, we can see that broadly the same individuals are perceived as sounding low and high SES. The main change is that speakers 29 and 37 are not perceived as sounding as broad in the social subset, although they are are still situated closer to the leaders than to the other laggers. We can see that the same speakers comprise the perceptually high SES groups in both spaces (22, 6, 28, 14, 18, 11, 35), and almost the same speakers comprise the perceptually low SES groups in both spaces (30, 2, 16,17, 31, 32, with 8 becoming more perceptually broad in the social subset). Speaker 13 is an outlier slow leader in both spaces, and the same fast leaders are grouped with perceptually low SES (20, 27, 7, 26) and perceptually high SES speakers (21,19,36,1,12) in both spaces. As such, the social differentiation between speakers is effectively captured by the full data MDS.
Click here to view code.
visualisation_df %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat, linetype=SES_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2,linetype=SES_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 == "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
alpha=0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
) +
geom_text(aes(label=SpeakerID), alpha=0.5)+
scale_fill_manual(values = c(c("#a541f7","#ffd117"))) +
scale_colour_manual(values = c(c("#a541f7", "#F9A742", "#ffd117","red"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "LL/Speed",
shape = "LL/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group",
title = "Full data"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "none", #change legend title font size
legend.text = element_text(size = 6))+
regression_df_social %>%
mutate(
PC1_tree2 = case_when(LeaderLaggerScore < -0.32 ~ 'Lagger', T ~ "Leader"),
speed_tree = case_when(ArticRate < (-0.2) ~ "Slow", T ~ "Fast"),
PC1_speed_tree = case_when(
PC1_tree2 == "Lagger" ~ "Lagger",
T ~ interaction(PC1_tree2, speed_tree)
),
) %>%
ggplot(aes(x = Rotated_D1_FC_social, y = Rotated_D2_FC_social)) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat, linetype=SES_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2, linetype=SES_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 == "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
alpha=0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
) +
geom_text(aes(label=SpeakerID), alpha=0.5)+
scale_fill_manual(values = c(c("#a541f7","#ffd117"))) +
scale_colour_manual(values = c(c("#a541f7", "#F9A742", "#ffd117","red"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "LL/Speed",
shape = "LL/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group",
title = "Social Subset"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "right", #change legend title font size
legend.text = element_text(size = 6))
7.5.1.1 Confirming relationship between social D1 and PC1
The code chunks below confirm that, as with the full data set, the social Label PC1 scores themselves are also predicted by the same variables, with some minor adjustments to the cutoff values for the leader-lagger scores and speed.
Click here to view code.
fit_Dim_SES <-
tree_spec %>%
fit(SES_Dim ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df_social)
fit_Dim_SES %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#ffd117", "#F9A742", "#be4580", "#e66f70", "#a541f7"),
cex = 1,
xflip = T
)
Figure 34 shows the overall interpretation of leaders and laggers remains stable for the social data set regardless of whether we use the cutoffs from the trees predicting D1 or PC1 for the social subset: the perceptually high SES speaker group consists predominantly of the same laggers, and the perceptually low SES group consists entirely of the same (mostly slow) leaders.
Click here to view code.
regression_df_social %>%
mutate(
PC1_tree = case_when(LeaderLaggerScore < -0.32 ~ 'Lagger', T ~ "Leader"),
speed_tree = case_when(ArticRate < (-0.2) ~ "Slow", T ~ "Fast"),
PC1_speed_tree = case_when(
PC1_tree == "Lagger" ~ "Lagger",
T ~ interaction(PC1_tree, speed_tree)
),
) %>%
ggplot(aes(x = Rotated_D1_FC_social, y = Rotated_D2_FC_social)) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat, linetype = SES_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2, linetype = SES_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 == "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
alpha = 0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
) +
geom_text(aes(label = SpeakerID), alpha = 0.5) +
scale_fill_manual(values = c(c("#a541f7", "#ffd117"))) +
scale_colour_manual(values = c(c("#a541f7", "#F9A742", "#ffd117", "red"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "LL/Speed",
shape = "LL/Speed",
fill = "Perceptual Group",
title = "D1 score cutoffs"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "none", #change legend title font size
legend.text = element_text(size = 6)) +
regression_df_social %>%
mutate(
PC1_tree = case_when(LeaderLaggerScore < -0.22 ~ 'Lagger', T ~ "Leader"),
speed_tree3 = case_when(ArticRate < (-0.18) ~ "Slow", T ~ "Fast"),
PC1_speed_tree = case_when(
PC1_tree == "Lagger" ~ "Lagger",
T ~ interaction(PC1_tree, speed_tree3)
),
) %>%
ggplot(aes(x = Rotated_D1_FC_social, y = Rotated_D2_FC_social)) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat != "Middle"),
aes(fill = SES_Dim_Cat, linetype = SES_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(fill = SES_Dim_Cat2, linetype = SES_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 == "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
alpha = 0.15
) +
geom_point(
data = . %>% filter(SES_Dim_Cat2 != "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
) +
geom_text(aes(label = SpeakerID), alpha = 0.5) +
scale_fill_manual(values = c(c("#a541f7", "#ffd117"))) +
scale_colour_manual(values = c(c("#a541f7", "#F9A742", "#ffd117", "red"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "LL/Speed",
shape = "LL/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group",
title = "Label PC cutoffs"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "right", #change legend title font size
legend.text = element_text(size = 6))
This is further upheld in the results of a random forest procedure predicting PC1; the importance of the leader-lagger continuum comes out even more strongly for PC1 scores than it does for D1 scores. In other words, the predictive power of the leader-lagger continuum is strongest when predicting the use of labels related to perceived socioeconomic status and accent broadness (PC1), which is systematically related to a speaker’s positive along Dimension 1.
Click here to view code.
set.seed(2)
PC1_social_rf_fit <-
random_f %>%
fit(SES_Dim ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df_social)
PC1_social_rf_fit %>%
extract_fit_engine() %>%
vip() +
theme_bw(base_size = 12) +
labs(y = "Importance (Random Forest)")
7.5.2 Predicting Social D2 with a regression tree
Figure 36 depicts the results of the regression tree predicting social D2 fit in the code chunk below. In contrast to D1, D2 for the social subset is not stably predicted by speed and pitch as D2 was for the full data set. Instead, the Back Vowel Configuration is represented twice in the tree.
Click here to view code.
set.seed(9)
fit_D2_social <-
tree_spec %>%
fit(
Rotated_D2_FC_social ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
MeanPitch2,
data = regression_df_social
)
fit_D2_social %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#09c9c3", "#8BDAB2", "#F9A742", "#F87233", "#f0496a"),
cex = 1,
xflip = T
)
For the sake of comparison, Figure 37 maps the same cutoffs from onto the MDS spaces from the full data set (also shown in Figure 20 ) and the social subset. We can see that some of the slow speakers (8, 32) have moved more into the ‘Young’ perceptual space from the ‘slow’ perceptual space, while some of the fast/high speakers (20 and 26) have moved from the ‘Fast/high’ perceptual space into the ‘Old’ perceptual space (29 only remains in the ‘Young’ perceptual space but has similarly moved along D2). There are also some fast and low speakers who have moved from the ‘Fast/high’ perceptual space into the ‘old’ perceptual space (15, 37, 3). From this, it is clear that perceptual age is being informed by something other than (or in addition to) the combination of speed and pitch that differentiates between speakers acoustically.
Click here to view code.
visualisation_df %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat != "Middle"),
aes(fill = Speed_Dim_Cat, linetype = Speed_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(fill = Speed_Dim_Cat2, linetype = Speed_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(shape = pitch_speed_tree,
colour = pitch_speed_tree),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat2 == "Middle"),
aes(shape = pitch_speed_tree,
colour = pitch_speed_tree),
size = 6,
stroke = 1,
alpha = 0.35
) +
geom_text(aes(label = SpeakerID)) +
scale_fill_manual(values = c("#f0496a", "#09c9c3")) +
scale_colour_manual(values = c("#f0496a", "#F9A742", "#09c9c3", "yellow")) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "Pitch/Speed",
shape = "Pitch/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group"
) +
theme_bw() +
coord_fixed() +
guides(shape = "none",
colour = "none") +
theme(legend.position = "bottom",
legend.text = element_text(size = 6)) +
regression_df_social %>%
mutate(
speed_tree = case_when(ArticRate < (-0.11) ~ "Slow", T ~ "Fast"),
pitch_tree2 = case_when(MeanPitch2 < (0.15) ~ "Low", T ~ "High"),
pitch_speed_tree = case_when(
speed_tree == "Slow" ~ "Slow",
speed_tree != "Slow" &
pitch_tree2 == "High" ~ "FastHigh",
T ~ "FastLow"
)
) %>%
ggplot(aes(x = Rotated_D1_FC_social, y = Rotated_D2_FC_social)) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat != "Middle"),
aes(fill = Age_Dim_Cat, linetype = Age_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat2 != "Middle"),
aes(fill = Age_Dim_Cat2, linetype = Age_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(Age_Dim_Cat2 != "Middle"),
aes(shape = pitch_speed_tree,
colour = pitch_speed_tree),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(Age_Dim_Cat2 == "Middle"),
aes(shape = pitch_speed_tree,
colour = pitch_speed_tree),
size = 6,
stroke = 1,
alpha = 0.35
) +
geom_text(aes(label = SpeakerID)) +
scale_linetype_manual(values = c("dashed", "dotdash")) +
scale_fill_manual(values = c("#8586fc", "#ff8164")) +
scale_colour_manual(values = c("#f0496a", "#F9A742", "#09c9c3", "yellow")) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "Pitch/Speed",
shape = "Pitch/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "right",
legend.text = element_text(size = 6))
While it is possible that the back-vowel configuration (BVC) is contributing to perceived age, when we map the cutoffs from the regression tree in Figure 36 onto the MDS spaces from the full and social data sets in Figure 38, there is not a discernible pattern. Most of the speakers in the perceptually ‘old’ group have BVC scores in middle of its range, not markedly higher or lower scores. Moreover, there are almost as many speakers with positive as negative scores (again classified as scores 0.5 SD +/- the mean) in the perceptually young/friendly group. It is therefore difficult to meaningfully interpret the space based on speakers’ Back-Vowel Configuration scores.
Click here to view code.
visualisation_df %>%
mutate(
BVC_tree = case_when(
BackVowelScore < (-0.63) ~ "Negative" ,
BackVowelScore >= 0.74 ~ "Positive",
T ~ "Middle"
)
) %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat != "Middle"),
aes(fill = Speed_Dim_Cat, linetype = Speed_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(fill = Speed_Dim_Cat2, linetype = Speed_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(shape = BVC_tree,
colour = BVC_tree),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat2 == "Middle"),
aes(shape = BVC_tree,
colour = BVC_tree),
size = 6,
stroke = 1,
alpha = 0.35
) +
geom_text(aes(label = SpeakerID)) +
scale_fill_manual(values = c(c("#f0496a", "#09c9c3"))) +
scale_colour_manual(values = c(c("#8586fc", "#9eee95", "#ff8164", "yellow"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "BVC Score",
shape = "BVC Score",
fill = "Perceptual Group",
linetype = "Perceptual Group"
) +
theme_bw() +
coord_fixed() +
guides(shape = "none",
colour = "none") +
theme(legend.position = "bottom",
legend.text = element_text(size = 6)) +
regression_df_social %>%
mutate(
BVC_tree = case_when(
BackVowelScore < (-0.63) ~ "Negative" ,
BackVowelScore >= 0.74 ~ "Positive",
T ~ "Middle"
)
) %>%
ggplot(aes(x = Rotated_D1_FC_social, y = Rotated_D2_FC_social)) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat != "Middle"),
aes(fill = Age_Dim_Cat, linetype = Age_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat2 != "Middle"),
aes(fill = Age_Dim_Cat2, linetype = Age_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(Age_Dim_Cat2 != "Middle"),
aes(shape = BVC_tree,
colour = BVC_tree),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(Age_Dim_Cat2 == "Middle"),
aes(shape = BVC_tree,
colour = BVC_tree),
size = 6,
stroke = 1,
alpha = 0.35
) +
geom_text(aes(label = SpeakerID)) +
scale_linetype_manual(values = c("dashed", "dotdash")) +
scale_fill_manual(values = c(YoungFriendly = "#ff8164", Old = "#8586fc")) +
scale_colour_manual(values = c(c("#8586fc", "#9eee95", "#ff8164", "yellow"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "BVC Score",
shape = "BVC Score",
fill = "Perceptual Group",
linetype = "Perceptual Group"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "right",
legend.text = element_text(size = 6))
7.5.2.1 Relationship between social D2 and PC2, and D2 and PC2 from full data set
Moreover, when we predict speaker’s label PC2 scores from the social data PCA with the same regression tree procedure in Figure 39, the predictors are not stable. Instead, articulation rate emerges as the most important predictor, mediated by speakers’ leader-lagger scores.
When the cutoffs from Figure 39 are mapped onto the MDS spaces from the full and social data in Figure 40 there is a slightly more coherent picture, in that all the very slow speakers, except speaker 8, are concentrated in the perceptually old group. Most of the perceptually young/friendly group is comprised of fast laggers, and most fast leaders are concentrated in the perceptually old group.
Click here to view code.
regression_df %>%
mutate(
PC1_tree = case_when(LeaderLaggerScore < -0.036 ~ 'Lagger', T ~ "Leader"),
speed_tree2 = case_when(ArticRate < (-0.5) ~ "Slow", T ~ "Fast"),
PC1_speed_tree = case_when(
speed_tree2 == "Slow" ~ "Slow",
speed_tree2 != "Slow" &
PC1_tree == "Lagger" ~ "FastLagger",
T ~ "FastLeader"
),
) %>%
ggplot(aes(x = Rotated_D1_FC, y = Rotated_D2_FC)) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat != "Middle"),
aes(fill = Speed_Dim_Cat, linetype = Speed_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(fill = Speed_Dim_Cat2, linetype = Speed_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat2 != "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(Speed_Dim_Cat2 == "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
alpha = 0.35
) +
geom_text(aes(label = SpeakerID)) +
scale_fill_manual(values = c(c("#f0496a", "#09c9c3"))) +
scale_colour_manual(values = c(c("#a541f7", "#F9A742", "#09c9c3", "yellow"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "Pitch/Speed",
shape = "Pitch/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group"
) +
theme_bw() +
coord_fixed() +
guides(shape = "none",
colour = "none") +
theme(legend.position = "bottom",
legend.text = element_text(size = 6)) +
regression_df_social %>%
mutate(
PC1_tree = case_when(LeaderLaggerScore < -0.036 ~ 'Lagger', T ~ "Leader"),
speed_tree2 = case_when(ArticRate < (-0.5) ~ "Slow", T ~ "Fast"),
PC1_speed_tree = case_when(
speed_tree2 == "Slow" ~ "Slow",
speed_tree2 != "Slow" &
PC1_tree == "Lagger" ~ "FastLagger",
T ~ "FastLeader"
),
) %>%
ggplot(aes(x = Rotated_D1_FC_social, y = Rotated_D2_FC_social)) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat != "Middle"),
aes(fill = Age_Dim_Cat, linetype = Age_Dim_Cat),
alpha = 0.35
) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat2 != "Middle"),
aes(fill = Age_Dim_Cat2, linetype = Age_Dim_Cat2),
alpha = 0.15
) +
geom_point(
data = . %>% filter(Age_Dim_Cat2 != "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(Age_Dim_Cat2 == "Middle"),
aes(shape = PC1_speed_tree,
colour = PC1_speed_tree),
size = 6,
stroke = 1,
alpha = 0.35
) +
geom_text(aes(label = SpeakerID)) +
scale_linetype_manual(values = c("dashed", "dotdash")) +
scale_fill_manual(values = c(YoungFriendly = "#ff8164", Old = "#8586fc")) +
scale_colour_manual(values = c(c("#a541f7", "#F9A742", "#09c9c3", "yellow"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "Pitch/Speed",
shape = "Pitch/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "right",
legend.text = element_text(size = 6))
It is, therefore, possible that leaders in change are perceived to sound older than fast laggers in change (along with very slow speakers). Counter-intuitively, the association of (slow) leaders with rurality may also carry connotations of social (rather than linguistic) conservatism. If laggers are less salient than leaders, it could be that rather than sounding explicitly ‘young’, such speakers do not carry association with conservatism and are therefore perceived as younger than leaders. More recent analyses of New Zealand English data also point to some reversals in the realisations of diphthongs (e.g., Ross, Ballard, and Watson 2022), meaning our ‘leaders’ may indeed sound older to our comparatively young listener sample.
7.5.2.2 Alternate interpretations of social D2
While it is possible that leaders and laggers are associated with different ages in addition to socioeconomic status, it is important to note that this association is not very strong, or stable within the MDS space (otherwise we would expect a stronger overlap between PC2 and D2, as observed for both the full data set, and PC1/D1 in both the full and social data). It is also worth considering the possibility that perceived age is associated with other independent variables not considered in the main analysis. For example, if we incorporate a measure of pause usage in the stimuli in the regression tree predicting D2 in the code chunk below, it emerges instead of the Back Vowel Configuration in Figure 41 (although the regression tree predicting PC2 remains unchanged).
Click here to view code.
set.seed(9)
fit_D2_social2 <-
tree_spec %>%
fit(
Rotated_D2_FC_social ~
BackVowelScore +
LeaderLaggerScore +
ArticRate +
PauseProp +
MeanPitch2,
data = regression_df_social
)
fit_D2_social2 %>%
extract_fit_engine() %>%
rpart.plot(
roundint = FALSE,
extra = 101,
box.palette = c("#09c9c3", "#8BDAB2", "#F9A742", "#F87233", "#f0496a"),
cex = 1,
xflip = T
)
One possible interpretation of the second dimension in the social subset is shown in Figure 42; the perceptually old group is comprised primarily of speakers who are either very slow or have a high use of pauses. Speakers who have a low use of pauses are concentrated in the perceptually young/friendly group.
Click here to view code.
regression_df_social_vis %>%
ggplot(aes(x = Rotated_D1_FC_social, y = Rotated_D2_FC_social)) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat == "YoungFriendly"),
aes(fill = Age_Dim_Cat, linetype = Age_Dim_Cat),
alpha = 0.15
) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat2 == "YoungFriendly"),
aes(fill = Age_Dim_Cat2, linetype = Age_Dim_Cat2),
alpha = 0.15
) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat2 == "Old"),
aes(fill = Age_Dim_Cat2, linetype = Age_Dim_Cat2),
alpha = 0.15
) +
geom_encircle(
data = . %>% filter(Age_Dim_Cat == "Old"),
aes(fill = Age_Dim_Cat, linetype = Age_Dim_Cat),
alpha = 0.15
) +
geom_point(
data = . %>% filter(Age_Dim_Cat2 != "Middle"),
aes(shape = pause_prop_time_tree,
colour = pause_prop_time_tree),
size = 6,
stroke = 1,
) +
geom_point(
data = . %>% filter(Age_Dim_Cat2 == "Middle"),
aes(shape = pause_prop_time_tree,
colour = pause_prop_time_tree),
size = 6,
stroke = 1,
alpha = 0.35
) +
geom_text(aes(label = SpeakerID)) +
scale_fill_manual(values = c(YoungFriendly = "#ff8164", Old = "#8586fc")) +
scale_colour_manual(values = c(c("#08c370", "#f94b5e", "#09c9c3", "yellow"))) +
labs(
x = "Rotated Dimension 1 (D1)",
y = "Rotated Dimension 2 (D2)",
colour = "Pitch/Speed",
shape = "Pitch/Speed",
fill = "Perceptual Group",
linetype = "Perceptual Group"
) +
theme_bw() +
coord_fixed() +
theme(legend.position = "right",
legend.text = element_text(size = 6))
The combined evidence points to a stable association between the leader-lagger continuum and perceived socioeconomic status, mediated by articulation rate. There is also a stable relationship between perceived and actual speed and pitch, with fast and pitch speakers particularly perceptually distinctive. The available evidence also indicates that the differentiation between fast/slow and low/high pitch speakers is driven primarily by the acoustic features listeners heard, rather than a strong association between these measures and age (within a single age category). While we cannot entirely rule out an association between speed, the leader-lagger vowels and perceived age, other factors are likely contributing to this perceptual distinction.
8 Packages used
We used R version 4.3.3 (R Core Team 2024) and the following R packages: colorspace v. 2.1.1 (Zeileis, Hornik, and Murrell 2009; Stauffer et al. 2009; Zeileis et al. 2020), corrplot v. 0.95 (Wei and Simko 2024), cowplot v. 1.1.3 (Wilke 2024), dials v. 1.4.0 (Kuhn and Frick 2025), e1071 v. 1.7.16 (Meyer et al. 2024), factoextra v. 1.0.7 (Kassambara and Mundt 2020), ggalt v. 0.4.0 (Rudis, Bolker, and Schulz 2017), ggcorrplot v. 0.1.4.1 (Kassambara 2023a), ggpubr v. 0.6.0 (Kassambara 2023b), gratia v. 0.10.0 (Simpson 2024), gt v. 1.0.0 (Iannone et al. 2025), here v. 1.0.1 (Müller 2020), infer v. 1.0.8 (Couch et al. 2021), kableExtra v. 1.4.0 (Zhu 2024), knitr v. 1.50 (Xie 2014, 2015, 2025), lattice v. 0.22.7 (Sarkar 2008), magick v. 2.8.6 (Ooms 2025), modeldata v. 1.4.0 (Kuhn 2024), nzilbb.vowels v. 0.4.1 (Wilson Black et al. 2023b), parsnip v. 1.3.1 (Kuhn and Vaughan 2025), patchwork v. 1.3.0 (Pedersen 2024), permute v. 0.9.7 (Simpson 2022), plotrix v. 3.8.4 (J 2006), randomForest v. 4.7.1.2 (Liaw and Wiener 2002), recipes v. 1.3.1 (Kuhn, Wickham, and Hvitfeldt 2025), rpart v. 4.1.24 (Therneau and Atkinson 2025), rpart.plot v. 3.1.2 (Milborrow 2024), rsample v. 1.3.0 (Frick et al. 2025), scales v. 1.4.0 (Wickham, Pedersen, and Seidel 2025), smacof v. 2.1.7 (de Leeuw and Mair 2009; Mair, Groenen, and de Leeuw 2022), tidymodels v. 1.3.0 (Kuhn and Wickham 2020), tidytext v. 0.4.2 (Silge and Robinson 2016), tidyverse v. 2.0.0 (Wickham et al. 2019), tune v. 1.3.0 (Kuhn 2025), vegan v. 2.6.10 (Oksanen et al. 2025), vip v. 0.4.1 (Greenwell and Boehmke 2020), workflows v. 1.2.0 (Vaughan and Couch 2025), workflowsets v. 1.1.0 (Kuhn and Couch 2024), yardstick v. 1.3.2 (Kuhn, Vaughan, and Hvitfeldt 2025).
References
Footnotes
Our preregistration specifies a ‘non-metric’ MDS. This typically means an ordinal MDS, but M-splines are also non-metric.↩︎
Note that the confidence bands for
NZProperthroughMediumPitchare very wide, indicating that the association of these variables with D2 is unlikely to hold for a different sample.↩︎The measure of creak proportion is from participants’ whole QB monologue and is therefore an indicator of their overall creakiness, not a measure of creak presence in the stimuli.↩︎
7.4 Social Label PCA
To see which label categories increase or decrease as the rotated social D1/D2 scores increase or decrease, we again implemented Principal Component Analysis (PCA). In this PCA the input data were D1/D2 scores and the proportion each social Level 1 label category (the most specific) was used to describe each stimuli (categories used by <10% of total participants were excluded).
Click here to view code.
Click here to view code.
7.4.1 Labels and Dimensions loaded onto PC1
Figure 29 depicts the variables loaded onto PC1 from the social label PCA. The results point to a similar pattern as for the full data for Dimension 1, where labels related to perceived socioeconomic status and accent strength are most strongly loaded onto PC1 (with social Dimension 1 scores). Specifically, the labels NZ Rural (accent), Strong NZ (accent), Low, Middle and High Socioeconomic status (SES) are the most strongly loaded speech-related labels, with NZ proper and British also marginally meeting our 90% threshold.
The loadings indicate that, as stimuli D1 scores increase, use of High and Middle SES, NZ proper and British labels decreases and use of Low SES, strong NZ (accent) and rural labels increases. Dimension 1 from the social MDS therefore appears to also primarily capture perceived social class and accent broadness.
Click here to view code.
7.4.2 Labels and Dimensions loaded onto PC2
Figure 30 depicts the variables loaded onto PC2 from the social label PCA. The results from PC2 point to a clearer split between socioeconomic labels and labels related to age for the social subset than for the full dataset. Specifically, labels related to both Young and Old age, are loaded onto label PC2, not just labels related to older age (as in Figure 10 ). Friendly and NZ Average also meet our 90% threshold. The social label PCA points to speakers labelled as younger as being unlikely to also be labelled as older, but perceived age is not systematically related to perceived accent strength and socioeconomic status. This points to the loading of ‘Young’ labels onto PC1 from the full data label PCA as not being very stable.
Click here to view code.
We extract speaker scores form the social PCA.
Click here to view code.