library(tidyverse)# tidyverse_2.0.0 library(lme4) # lme4_1.1-37 library(lmerTest) # lmerTest_3.1-3 library(emmeans) # emmeans_1.11.1 library(patchwork) # patchwork_1.3.2 library(performance) # performance_0.17.0 source("~/Library/RScripts/recontrast.R") # https://www.ling.upenn.edu/~joseff/scripts/recontrast.R #################################### ## Read in data #################################### data <- read_csv("Campbell-Kibler_Clopper_data.csv") %>% mutate(Phone = recontrast(factor(Phone)), Condition = relevel(factor(Condition), ref = "NoInstruction"), SoundFile = factor(SoundFile)) stimuli <- read_csv("Campbell-Kibler_Clopper_stimulus_data.csv") #################################### ## Main model #################################### modelDiff.lmer <- data %>% lmer(f = FormantDiff ~ Phone * Block * Condition + (1 + Phone + Block | SoundFile) + (1 + Block | Word), control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(modelDiff.lmer) anova(modelDiff.lmer) r2_nakagawa(modelDiff.lmer) modelDiff.emm <- emmeans(modelDiff.lmer, pairwise ~ Block | Condition * Phone, lmer.df="satterthwaite", lmerTest.limit = 11000) #################################### ## Main model in Bark (not presented in paper) #################################### modelDiffBark.lmer <- data %>% mutate(LowerFBark = (26.81 * LowerF)/(1960 + LowerF) - 0.53, HigherFBark = (26.81 * HigherF)/(1960 + HigherF) - 0.53, FormantDiffBark = HigherFBark - LowerFBark) %>% lmer(f = FormantDiffBark ~ Phone * Block * Condition + (1 + Phone + Block | SoundFile) + (1 + Block | Word), control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(modelDiffBark.lmer) anova(modelDiffBark.lmer) modelDiffBark.emm <- emmeans(modelDiffBark.lmer, pairwise ~ Block | Condition * Phone, lmer.df="satterthwaite", lmerTest.limit = 11000) #################################### ## Gender model #################################### modelDiffGender.lmer <- data %>% filter(Gender == "F" | Gender == "M") %>% lmer(f = FormantDiff ~ Phone * Block * Condition * Gender + (1 + Phone + Block | SoundFile) + (1 + Block + Gender | Word), control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(modelDiffGender.lmer) anova(modelDiffGender.lmer) modelDiffGender.emm <- emmeans(modelDiffGender.lmer, pairwise ~ Block | Gender * Phone * Condition, lmer.df="satterthwaite", lmerTest.limit=11000) #################################### ## Social salience models #################################### # We flip the sign on BATH measurements here to align convergence effects across the two phones # In other words bigger = more convergence in this model, regardless of phone commentary.lmer <- data %>% filter(Phone == "bath" | Phone == "near") %>% mutate(FormantDiff = ifelse(Phone == "bath", -1*FormantDiff, FormantDiff)) %>% lmer(f = FormantDiff ~ Block * Condition * NoticePhone + Phone * Condition * NoticePhone + Phone * Block * Condition + (1 + Phone * Block | SoundFile) + (1 + Block | Word), control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(commentary.lmer) anova(commentary.lmer) commentary.emm <- emmeans(commentary.lmer, pairwise ~ Block | Condition* NoticePhone, lmer.df="satterthwaite", lmerTest.limit = 11000) # TRAP time effect, in Imitate condition only modelTRAP_timeEffect.lmer <- data %>% filter(Phone == "trap", Condition == "Imitate") %>% mutate(TrialTime = WordStart - 200*(Block == "shadow")) %>% lmer(f = FormantDiff ~ Block*TrialTime+ (1 + Block+scale(TrialTime) | SoundFile) + (1+Block|Word), control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(modelTRAP_timeEffect.lmer) anova(modelTRAP_timeEffect.lmer) modelTRAP_timeEffect.emm <- emmeans(modelTRAP_timeEffect.lmer, ~ TrialTime | Block, lmer.df="satterthwaite", lmerTest.limit = 11000) #################################### ## Figure 1 #################################### theme_paper <- theme_gray() + theme(text = element_text (size = 15), panel.background = element_rect(fill = "white"), panel.grid.major = element_blank()) + theme(plot.title = element_text(hjust = 0.5)) theme_set(theme_paper) cexFactor <- 8 talkerPlot <- stimuli %>% filter(Phone != "near") %>% group_by(Phone) %>% summarise(F2 = mean(HigherF), F1 = mean(LowerF)) %>% ggplot() + geom_text(aes(x = F2, y = F1, label = str_to_upper(Phone)), show.legend = FALSE, color = "black", cex = cexFactor) + ylim(1100, 445) + xlim(2700, 1100) womensMeans <- data %>% filter(Gender == "F", Block == "baseline", Phone != "near") %>% mutate(Phone = fct_collapse(Phone, TRAP = c("trap", "bath"))) %>% group_by(Phone) %>% summarise(F2 = mean(HigherF), F1 = mean(LowerF)) png("Figure1.png") talkerPlot + geom_text(data = womensMeans, aes(x = F2, y = F1, label = str_to_upper(Phone)), show.legend = FALSE, color = "darkgrey", cex = cexFactor) + labs(x = "F2 (Hz)", y = "F1 (Hz)") + theme(panel.grid.major = element_line(color = "gray"), axis.ticks.x = element_blank(), axis.ticks.y =element_blank()) dev.off() #################################### ## Figure 2 #################################### subjMeans <- data %>% group_by(SoundFile, Phone, Condition, Block) %>% summarise(subjMeanDiff = mean(FormantDiff)) %>% ungroup() %>% mutate(ConditionBlock = paste(Condition, Block), ConditionBlock = factor(ConditionBlock, levels = c("AvoidImitate baseline", "AvoidImitate shadow", "NoInstruction baseline", "NoInstruction shadow", "Imitate baseline", "Imitate shadow"))) # This controls the gap between the line and the text label for the model talker ##### Bath plot bathPlot <- subjMeans %>% filter(Phone == "bath") %>% group_by(Condition, Block, ConditionBlock) %>% summarise(meanDiff = mean(subjMeanDiff), sdDiff = sd(subjMeanDiff), n = n()) %>% mutate(seDiff = sdDiff/sqrt(n)) %>% ggplot(aes(x = ConditionBlock, color = Block)) + geom_point(aes(y = subjMeanDiff), data = subset(subjMeans, Phone == "bath"), color = "lightgray", cex = 0.75) + geom_point(aes(y = meanDiff), cex = 2) + geom_errorbar(aes(ymin = meanDiff - seDiff, ymax = meanDiff + seDiff)) + scale_color_manual(values = c("black", "darkgray"), guide = "none") + labs(title = "BATH", x = "Condition", y = "F2-F1 @ 50% (Hz)") + annotate(geom = "text", x = c(1.5, 3.5, 5.5), y = min(subset(subjMeans, Phone == "bath")$subjMeanDiff) - 60, label = c("Avoid", "No Instr.", "Imitate")) + theme(panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x=element_blank()) + geom_hline(yintercept = 495, lty = 2) + annotate(geom = "text", x = 2, y = 555, label = "Model talker") ##### Trap plot trapPlot <- subjMeans %>% filter(Phone == "trap") %>% group_by(Condition, Block, ConditionBlock) %>% summarise(meanDiff = mean(subjMeanDiff), sdDiff = sd(subjMeanDiff), n = n()) %>% mutate(seDiff = sdDiff/sqrt(n)) %>% ggplot(aes(x = ConditionBlock, color = Block)) + geom_point(aes(y = subjMeanDiff), data = subset(subjMeans, Phone == "trap"), color = "lightgray", cex = 0.75) + geom_point(aes(y = meanDiff), cex = 2) + geom_errorbar(aes(ymin = meanDiff - seDiff, ymax = meanDiff + seDiff)) + scale_color_manual(values = c("black", "darkgray"), guide = "none") + labs(title = "TRAP", x = "Condition", y = "F2-F1 @ 50% (Hz)") + annotate(geom = "text", x = c(1.5, 3.5, 5.5), y = min(subset(subjMeans, Phone == "trap")$subjMeanDiff) - 60, label = c("Avoid", "No Instr.", "Imitate")) + theme(panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x=element_blank()) + geom_hline(yintercept = 1044, lty = 2) + annotate(geom = "text", x = 2, y = 1104, label = "Model talker") ##### Kit plot kitPlot <- subjMeans %>% filter(Phone == "kit") %>% group_by(Condition, Block, ConditionBlock) %>% summarise(meanDiff = mean(subjMeanDiff), sdDiff = sd(subjMeanDiff), n = n()) %>% mutate(seDiff = sdDiff/sqrt(n)) %>% ggplot(aes(x = ConditionBlock, color = Block)) + geom_point(aes(y = subjMeanDiff), data = subset(subjMeans, Phone == "kit"), color = "lightgray", cex = 0.75) + geom_point(aes(y = meanDiff), cex = 2) + geom_errorbar(aes(ymin = meanDiff - seDiff, ymax = meanDiff + seDiff)) + scale_color_manual(values = c("black", "darkgray"), guide = "none") + labs(title = "KIT", x = "Condition", y = "F2-F1 @ 50% (Hz)") + annotate(geom = "text", x = c(1.5, 3.5, 5.5), y = min(subset(subjMeans, Phone == "kit")$subjMeanDiff) - 60, label = c("Avoid", "No Instr.", "Imitate")) + theme(panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x=element_blank()) + geom_hline(yintercept = 2168, lty = 2) + annotate(geom = "text", x = 2, y = 2235, label = "Model talker") ##### Price plot pricePlot <- subjMeans %>% filter(Phone == "price") %>% group_by(Condition, Block, ConditionBlock) %>% summarise(meanDiff = mean(subjMeanDiff), sdDiff = sd(subjMeanDiff), n = n()) %>% mutate(seDiff = sdDiff/sqrt(n)) %>% ggplot(aes(x = ConditionBlock, color = Block)) + geom_point(aes(y = subjMeanDiff), data = subset(subjMeans, Phone == "price"), color = "lightgray", cex = 0.75) + geom_point(aes(y = meanDiff), cex = 2) + geom_errorbar(aes(ymin = meanDiff - seDiff, ymax = meanDiff + seDiff)) + scale_color_manual(values = c("black", "darkgray"), guide = "none") + labs(title = "PRICE", x = "Condition", y = "F2-F1 @ 35% (Hz)") + annotate(geom = "text", x = c(1.5, 3.5, 5.5), y = min(subset(subjMeans, Phone == "price")$subjMeanDiff) - 60, label = c("Avoid", "No Instr.", "Imitate")) + theme(panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x=element_blank()) + geom_hline(yintercept = 428, lty = 2) + annotate(geom = "text", x = 2, y = 488, label = "Model talker") ##### Near plot nearPlot <- subjMeans %>% filter(Phone == "near") %>% group_by(Condition, Block, ConditionBlock) %>% summarise(meanDiff = mean(subjMeanDiff), sdDiff = sd(subjMeanDiff), n = n()) %>% mutate(seDiff = sdDiff/sqrt(n)) %>% ggplot(aes(x = ConditionBlock, color = Block)) + geom_point(aes(y = subjMeanDiff), data = subset(subjMeans, Phone == "near"), color = "lightgray", cex = 0.75) + geom_point(aes(y = meanDiff), cex = 2) + geom_errorbar(aes(ymin = meanDiff - seDiff, ymax = meanDiff + seDiff)) + scale_color_manual(values = c("black", "darkgray")) + labs(title = "NEAR", x = "Condition", y = "F3-F2 @ 80% (Hz)") + annotate(geom = "text", x = c(1.5, 3.5, 5.5), y = min(subset(subjMeans, Phone == "near")$subjMeanDiff) - 60, label = c("Avoid", "No Instr.", "Imitate")) + theme(panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x=element_blank()) + geom_hline(yintercept = 1019, lty = 2) + annotate(geom = "text", x = 2, y = 1100, label = "Model talker") ## Put them all together png("Figure2.png", width = 960) design <- "AABBCC #DDEE#" bathPlot + trapPlot + kitPlot + pricePlot + nearPlot + plot_layout(design = design) dev.off()