Ongoing exposure to an ambient language continues to build implicit knowledge across the lifespan
June 10, 2024
1 Introduction
This R Markdown file contains supplementary information and all the code used for data exclusion, analysis, and plotting. This document is split into five sections. Section 2 describes the data preparation process, including the dataset structure, the factors associated with the stimuli, and how they were calculated. It also details the experimental methodology and participant demographics. Section 3 describes and provides the code for the analysis of the Adolescent Data for both the Wellformedness Rating Task and the Word Identification Task. Section 4 describes and supplies the code for the analysis of all datasets. It also details the comparison of adult and adolescent data to examine the role of continued exposure into adulthood in increasing the size of the proto-lexicon. Section 5 describes and provides the code for the analysis of the Adult Data in regard to the role of exposure on the proto-lexicon.
library(pacman)
p_load(knitr, tidyverse, ordinal, effects, MASS, kableExtra, egg, lme4, ggpubr, emmeans,car)
#library(knitr)
#library(tidyverse)
#library(ordinal)
#library(effects)
#library(MASS)
#library(kableExtra)
#ibrary(egg)
opts_chunk$set(echo=TRUE, message=FALSE, warning=FALSE, fig.show='hold', results='hold')
# Local functions ----
# A function to calculate confidence intervals
CI <- function (x, ci = 0.95)
{
a <- mean(x)
s <- sd(x)
n <- length(x)
error <- qt(ci + (1 - ci)/2, df = n - 1) * s/sqrt(n)
return(c(upper = a + error, mean = a, lower = a - error))
}
# A function to center variables
c. <- function (x) scale(x, scale = FALSE)
# Custom function to format ordinal regression results in a nice table
clm_table <- function(mod, digits=3, ...) {
table <- summary(mod)$coefficients %>%
data.frame() %>%
rlang::set_names(c("beta", "se", "z", "p")) %>%
rownames_to_column("parameter") %>%
mutate(
parameter = parameter %>%
str_replace_all(
.,
str_c("(", str_c(names(attr(mod$terms, "dataClasses"))[attr(mod$terms, "dataClasses") %in% c("factor", "logical", "character")], collapse="|"), ")"),
"\\1 = "
) %>%
str_replace_all(
.,
"(?<=^|\\:)(?:c\\.\\()([^:]+)(?:\\))(?=$|\\:)",
"\\1 (centered)"
) %>%
str_replace_all(., fixed(":"), " × "),
significance = case_when(
p < 0.001 ~ "***",
p < 0.01 ~ "**",
p < 0.05 ~ "\\*",
p < 0.1 ~ ".",
TRUE ~ ""
),
is_threshold = parameter %in% names(mod$alpha)
) %>%
mutate_at(c("z", "p", "significance"), ~ ifelse(is_threshold, NA, .)) %>%
arrange(is_threshold) %>%
mutate(
type = c("Effects", rep("", n()-sum(is_threshold)-1), "Thresholds", rep("", sum(is_threshold)-1))
) %>%
dplyr::select(type, parameter, beta, se, z, p, significance) %>%
mutate(
p = ifelse(p<0.001, "<0.001", format(round(p, 3), nsmall=3))
) %>%
kable(digits=3, escape=F, col.names=c("", "Parameter", "Estimate", "Std. Error", "$z$", "$p$", ""), align="llrrrrl", ...) %>%
column_spec(1, bold=TRUE) %>%
kable_styling()
return(table)
}
options(knitr.kable.NA = '')
# Custom function to compute summary data for an ordinal regression model, to be plotted. The "type" argument specifies whether the data are for a latent variable plot, a distributional plot, or a mean rating plot.
clm_plotdat = function(mod, focal.predictors, type="latent", response.coding=NULL, ...) {
# For latent variable plot, draw straight from Effect
if (type=="latent") {
latentdat <-
Effect(focal.predictors, mod, latent=TRUE, ...) %>%
as.data.frame() %>%
dplyr::select(-se) %>%
rename("pred"=fit, "lci"=lower, "uci"=upper)
return(latentdat)
}
# Get distributional data (used for both of other plots)
distdat <-
Effect(focal.predictors, mod, ...) %>%
as_tibble() %>%
dplyr::select(all_of(focal.predictors), matches("^(?:(?:L|U)\\.)?prob")) %>%
pivot_longer(
cols = -all_of(focal.predictors),
names_to = c(".value", "response"),
names_pattern = "^((?:(?:L|U)\\.)?prob)\\.(.+)$"
) %>%
rename("pred"=prob, "lci"=L.prob, "uci"=U.prob) %>%
mutate(
response = rep_len(mod$y.levels, n())
)
if (type=="dist") return(distdat)
if (type=="mean") {
# Get threshold values
thresholds = c(-Inf, as.numeric(mod$alpha), Inf)
# Get response coding
if (is.null(response.coding)) {
response.coding = 1:(length(mod$y.levels))
names(response.coding) = mod$y.levels
}
# Build latent variable data based on distributional data
# Note: this is a HACK since drawing latent variable data from Effect doesn't account for threshold variability. The confidence intervals assume fixed thresholds, which obscures patterns in variation of the thresholds; CIs are thus not correct (but are better than what would be obtained from direct transformation of the latent variable from Effect)
# Get distribution-based latent variable data, assuming fixed thresholds
latentdat <- distdat %>%
group_by_at(focal.predictors) %>%
summarise(
pred = qlogis(1-pred[1]) + thresholds[2],
lci = qlogis(1-uci[1]) + thresholds[2],
uci = qlogis(uci[length(uci)]) + thresholds[length(thresholds)-1]
) %>%
as.data.frame()
# Get probabilities for each response option
meandat <- latentdat
for (i in 1:length(mod$y.levels)) {
meandat[, paste0("(", mod$y.levels[i], ").prob")] = plogis(latentdat[, "pred"] - thresholds[i]) - plogis(latentdat[, "pred"] - thresholds[i+1])
meandat[, paste0("(", mod$y.levels[i], ").lci")] = plogis(latentdat[, "lci"] - thresholds[i]) - plogis(latentdat[, "lci"] - thresholds[i+1])
meandat[, paste0("(", mod$y.levels[i], ").uci")] = plogis(latentdat[, "uci"] - thresholds[i]) - plogis(latentdat[, "uci"] - thresholds[i+1])
}
meandat <- meandat %>%
dplyr::select(-pred, -uci, -lci) %>%
pivot_longer(
cols = -all_of(focal.predictors),
names_to = c("response", ".value"),
names_pattern = "^\\((.+)\\)\\.([^.]+)$"
) %>%
left_join(
data.frame(response=names(response.coding), code=response.coding, stringsAsFactors=FALSE),
by = "response"
) %>%
group_by_at(focal.predictors) %>%
summarise(
pred = sum(prob * code),
lci = sum(lci * code),
uci = sum(uci * code)
) %>%
ungroup()
return(meandat)
}
}
# Custom function to get the information about thresholds, in a nice format for plotting
clm_thresholds = function(mod) {
thresholds <- summary(mod)$coefficients %>%
as.data.frame() %>%
rownames_to_column() %>%
dplyr::select(1:3) %>%
set_names(c("label", "value", "se")) %>%
filter(label %in% names(mod$alpha)) %>%
mutate(
lci = value - 1.96*se,
uci = value + 1.96*se
) %>%
dplyr::select(-se)
return(thresholds)
}
# A function to add phonotactic scores from a file to the dataset
add_scores = function(data, scoreFile, scoreName) {
scored_data <- data %>%
left_join(
read_csv(scoreFile, col_types=cols_only(item=col_character(),
logprob=col_double())
)%>%
transmute(
item = item,
!!scoreName := logprob / (nchar(item) + 1)
)
, by="item"
)
return(scored_data)
}
# A function to run shell scripts on Windows or UNIX-based OS
run_shell = function(cmd) {
if (Sys.info()[['sysname']] == "Windows") {
return(shell(cmd, flag="", ignore.stdout=TRUE, ignore.stderr=TRUE))
} else {
return(system(cmd, ignore.stdout=TRUE, ignore.stderr=TRUE))
}
}
# A function to highlight the minimum value in a column, for presentation in a kable
# Highlighting is accomplished by coloring red
highlight_min = function(col, digits=1) {
rounded <- round(col, digits)
i <- which(rounded==min(rounded))
highlighted <- cell_spec(format(rounded, nsmall=digits), align="r")
highlighted[i] <- cell_spec(format(rounded[i], nsmall=digits), align="r", color="red")
return(highlighted)
}
# A function to display a table of numbers
display_table = function(dat, caption, digits=3, label=NA, highlight=c()) {
table <- dat %>%
mutate_at(highlight, ~ highlight_min(., digits=digits)) %>%
kable(caption=caption, digits=digits, escape=F, label=label) %>%
kable_styling()
return(table)
}
# A table for emmeans results
em_table = function(mod, ...){
table <- summary(mod) %>% mutate(across(where(is.numeric), round, 3)) %>% mutate(
significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "\\*",
p.value < 0.1 ~ ".",
TRUE ~ ""),
p.value = ifelse(p.value<0.001, "<0.001", format(round(p.value, 3), nsmall=3))) %>%
dplyr::select(-df)%>%
# kable(digits = 3,col.names=c("Parameter","Estimate","Std. Error","lower CL","upper CL","$z$", "$p$", ""),align="lrrrrrrl") %>% kable_styling()
#return(table)
return()
}2 Data Preparation and Method
2.1 Material
Adult Data (Data from Panther et al. 2023)
wellFormThe Well-Formedness Rating Task Data For the Wellformedness Rating Task, 40 non-word stimuli were randomly sampled from each of the high, medium, and low phonotactic score bins, for a total of 120 stimuli.wordIdentThe Word Identification Task Data
For the Identification Task, 20 stimulus pairs were randomly sampled from each of the high, mid, and low frequency bins. This yielded 60 real words and 60 matched non-words, for a total of 120 stimuli. Each participant in the experiment received a random sample, and each sample was shuffled into a random order.
See Panther et al. (2023) for detailed material and methods.
Adolescent Data (Our adolescent experiment)
wellForm20 non-word stimuli were sampled from each of the high, medium, and low phonotactic score bins, for a total of 60 stimuli.wordIdent60 real words and 60 matched non-words were sampled from the high and mid frequency bins, for a total of 120 stimuli. Each participant in the experiment received the same sample but the sample was shuffled into a random order for each participant.
The task was made shorter to make it more manageable, and that we did some post-hoc modelling of previous data to establish a length that would still give sufficient power. We used the same words for all adolescent to try and maximize comparability across participants, in case there were clear age effects.
2.2 Task Procedure
Wellformedness TaskParticipants were presented with each stimulus one at a time and rated it on a 5-point Likert scale based on how Māori-like they perceived it to be. A rating of 1 was labeled as “non-Māori-like non-word”, while a rating of 5 was labeled as “highly Māori-like non-word”.Identification TaskParticipants were presented with each stimulus one at a time and rated it on a 5-point Likert scale based on their confidence that it was a Māori word. A rating of 1 was labeled as “confident that this is NOT a Māori word”, while a rating of 5 was labeled as “confident that this is a Māori word”.
2.3 Data Filtering
Adult Data
In the adult dataset, we used answers in the post-experiment questionnaire to filter participants from each data; 226 participants completed the experiment (Identification and Wellformedness Tasks), of whom 43 participants were removed. We filtered out those who met any of the following criteria:
- not completing the post-experiment questionnaire
- indicating that they could speak or understand Māori at least “fairly well”
- indicating that they had studied Linguistics
- indicating that they spoke a Polynesian language
- completing the experiment too quickly, as indicated by having a median reaction time that was below two standard deviations of the grand mean of participant-wise median reaction times.
- indicating that they had lived outside NZ for longer than a year since they were 7
- indicating that their first language was not English
- indicating that they learned their first language outside New Zealand and have lived in their current location for a short period of time
Thus, the adult data consisted of 183 participants for the Identification Task. Firstly, we extracted two groups from this sample based on participants’ self-reported Māori Education: people who have never learnt Māori formally (n = 74) and people who have been exposed to Māori at a primary/intermediate school (n = 78). Three participants who were not in the Wellformedness Task Data were removed from Identification Task Data (n = 149 [never: n = 73, school: n = 76]).
# adult non-Māori-speaking NZers data
NZwellForm <- read_csv("docs/NZwellFormTask.csv")
#names(NZwellForm) #set Including Identification Task score
NZwellForm$workerId<- paste0("nz", NZwellForm$workerId)
#length(unique(NZwellForm$word)) #521
# the subset of the participants: select people do not learn Māori at all
NZwellFormNonMao<-NZwellForm %>% filter(maoriEducation %in% c("0","1"))
#length(unique(NZwellFormNonMao$workerId)) #149
#tableEdu<-NZwellFormNonMao%>% dplyr::select(workerId,maoriEducation) %>% distinct_all () %>% group_by(maoriEducation) %>% summarise(n()) # 0: 73, 1:76# read stimuli
stimuli<- read_csv("docs/identTaskStimuli.csv")
# how many stimuli?
# length(unique(NZident$word)) #544
# adult non-Māori-speaking NZers data
NZident <- read_csv("docs/NZidentTask.csv")
#length(unique(NZident$workerId)) 183
NZident$workerId<- paste0("nz", NZident$workerId)
# check how many participants didn't have any formal Māori education
NZident0Mao <- NZident %>% dplyr::select(workerId,maoriEducation) %>% distinct_all () %>% group_by(maoriEducation) %>% summarise(n())
# Q3: What is the highest level of education at which you have studied te reo Māori?", "school = 1"/ "high school = 2 "/ "undergraduate level = 3 "/ "At postgraduate level in university = 4"/"never =0"
# select people do not learn Māori at all
NZidentNonMao<-NZident %>% filter(maoriEducation %in% c("0","1"))
#length(unique(NZidentNonMao$workerId)) # 0:74,1:78 = 152
# remove workerId is not in the wellformedness task
NZidentNonMao<-NZidentNonMao[NZidentNonMao$workerId%in%NZwellFormNonMao$workerId,]
#length(unique(NZidentNonMao$workerId)) # 0:73,1:76 = 149Results in plots (Figure1) between two groups in raw data are not much different from each other. Therefore we use these two groups.
# check how different between people who have never learnt Māori formally and exposed at a primary school)
# Plot Ident Task with 149 participants: 73 pax(0 knowledge), 76 pax (1: learned at a school)
NZidentNonMao$word = as.character(NZidentNonMao$word)
Encoding(NZidentNonMao$word) = "UTF-8"
colnames(NZidentNonMao)[10] <- "freq.c"
MaoEducationPlot =NZidentNonMao %>%
group_by(maoriEducation, p.scr, type, freq.c) %>%summarise(
mean_response = mean(as.numeric(enteredResponse)),
n = n()
) %>%
ungroup() %>% mutate(freq.c=fct_relevel(freq.c,c("low", "mid","high"))) %>%
ggplot(., aes(x=p.scr, y=mean_response, color=type)) +
geom_point(aes(shape=type), alpha=0.3, size=4) +
geom_smooth(aes(fill=type), method="lm", formula="y~x", alpha=0.6, size=1) +
facet_grid(maoriEducation~ freq.c, labeller="label_both", switch="y") +
xlab("Phonotactic score") +
ylab("Mean rating (per stimulus)") +
scale_shape_manual(name="Stimulus type ",
values = c("nonword" = 1, "word" = 2),
guide = guide_legend(reverse = TRUE)) +
scale_color_manual(name="Stimulus type ",
values = c("nonword" = "black", "word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="Stimulus type ",
values = c("nonword" = "gray70", "word" = "dodgerblue1"),
guide = guide_legend(reverse = TRUE)) +
theme_bw() +
ggtitle("Effect of Māori Study") +
ylim(1, 5) +
theme(
panel.grid = element_blank(),
)
MaoEducationPlot
remove(MaoEducationPlot)
Figure 1:Group comparison between people who never studied Māori (n = 73) and people exposed at a primary school (n = 76). Mean wordhood confidence rating vs. phonotactic score for each stimulus for real words and non-words by frequency bin. Lines show correlations within each bin, for each stimulus type.
Adolescent Data
Step 1 The full dataset contains 58 participants. We
removed one participant that did not display much variation in their
responses, as identified by the standard deviation of their responses
(i.e. the 1-5 rating) being below two standard deviations of the grand
mean of participant-wise response standard deviations. Further we
removed two participants based on their median reaction time which was
below 2 standard deviations of the grand mean of participants-wise
median reaction time, indicating they finished the experiment too
quickly.
Step 2 For the next step, we removed participants based
on a demographic survey and a word-definition task. In a similar manner
as Adult Data, we extracted participants who have little knowledge of
the Māori language from Adolescent Data: Participants who provided a
definition for less than two words in a simple definition task (none:
n = 34, 1: n = 12). As a result, we have 46
adolescents.
Questions used in the demographic survey are:
What is your age?
What is your gender?
What is your ethnicity?
Have you ever been to school outside of New Zealand?
What language do you most often speak at home?
Can you speak or understand Māori?
- Yes, fluently
- Yes, fairly well
- Yes, but not very well (I can only understand or talk about simple things)
- No, only a few words or sentences
- Not at all
- Here is a list of 10 Māori words. For each word, if you know what it means, please write the meaning in English. If you don’t know what it means, just write ‘Don’t know’. Don’t worry if you don’t know the words – we’re not expecting you to! We just want to get some sense of how many words you might know.
- Word 1 pango (black)
- Word 2 ako (to learn)
- Word 3 oma (to run)
- Word 4 pō (night)
- Word 5 titiro (to look at)
- Word 6 pānui (to read)
- Word 7 kōtiro (girl)
- Word 8 rorohiko (computer)
- Word 9 hiamoe (sleepy)
- Word 10 upoko (head)
As a result, we have 195 participants from two groups for our statistical analysis.
# school kids data
SchoolIdent <- read_csv("docs/SchoolIdentTask.csv")
#names(SchoolIdent)
SchoolIdent$dataset<-"child"#(i.e.,"adolescent");length(unique(SchoolIdent$workerId)) #55
# check how many participants didn't have any formal Māori education
SchoolIdent0Mao<-SchoolIdent %>% dplyr::select(workerId, correct) %>% distinct_all () %>% group_by(correct) %>% summarise(n())
# School Ident data without vocabulary knowledge
SchoolIdent0M<-SchoolIdent %>% filter(correct %in% c("0","1"))
#length(unique(SchoolIdent0M$workerId)) # 46#length(unique(NZident$word)) #544
# Since two cohort are not much different, we use this subset data.
# file with 120items to match other data
NZident120<-NZidentNonMao[NZidentNonMao$word%in%stimuli$item,]
#length(unique(NZident120$word)) #120
# check how many pairs are used.
pair_counts<-NZident120%>% dplyr::select(workerId,pair.number,type) %>% distinct_all () %>% group_by(pair.number,type) %>% summarise(n())
#summary(pair_counts)#(minimum 8pairs, max57pairs)
#table(NZident120$workerId) #(minimum 20 stimuli, max 42 stimuli)
NZident120$dataset<-"adult"
#list(unique(NZident120$word))
#length(unique(NZident120$workerId)) #149
# stimulus lists
SchoolIdentList<-SchoolIdent %>% dplyr::select(word) %>% distinct(word) %>% arrange(word)
# original data (adult data)
OrigIdent<-NZident120 %>% dplyr::select (1:29, 60)
#names(OrigIdent)
setOriginalIdent <- OrigIdent %>% dplyr::select (word,poss.neighbors,nm.nbr) %>% distinct_all()
# school kids data
KidIdent<-SchoolIdent0M%>% dplyr::select (1,3:7, 9, 11:14,16:33,54) %>%
dplyr::rename(n.phon = n.phonemes, nm.nbr = norm.neighbors, p.scr = score.shortv, freq.c= freq.category)
#names(KidIdent)
remove(NZident, pair_counts,setOriginalIdent, SchoolIdentList)
#Comparing column names of two files
#fun <- function(OrigIdent,KidIdent) {
# if(all(names(OrigIdent) %in% names(KidIdent))) {
# cat('\nNo-Mismatch\n')
# KidIdent[names(OrigIdent)]
# }
# else {
# cat('\nNo cols : ', setdiff(names(OrigIdent), names(KidIdent)))
# cat('\nExtra cols : ', setdiff(names(KidIdent), names(OrigIdent)))
# }
#}
#fun(OrigIdent, KidIdent) #No-Mismatch
# bind three files except set 1 in order to extract Identification Task score for each participants
#length(unique(OrigIdent$workerId)) #149
#length(unique(KidIdent$workerId)) #46
identTask<-rbind(OrigIdent,KidIdent)
#total
#length(unique(identTask$workerId)) #195
#length(unique(identTask$word)) #120
remove(SchoolIdent0Mao, stimuli)# school kids data
SchoolWellForm <- read_csv("docs/SchoolWellFormTask.csv")
#names(SchoolWellForm)
SchoolWellForm$dataset<-"child" #(i.e.,adolescent)
# get participants match to Indent task
SchoolWellForm0M<-SchoolWellForm[SchoolWellForm$workerId %in% KidIdent$workerId,]
#length(unique(SchoolWellForm0M$workerId)) #46
nb_res <- aggregate(SchoolWellForm0M$trialNumber, by=list(SchoolWellForm0M$workerId), FUN=length)
#unique(nb_res$x);names(nb_res)<-c("workerId", "responses")
remove(KidIdent,OrigIdent)
# stimulus lists
SchoolList<-SchoolWellForm %>% dplyr::select(word) %>% distinct(word) %>% arrange(word)
# original non-Māori-speaking NZers data
#NZwellForm <- read_csv("NZwellFormTask.csv")
#names(NZwellForm) #set Including Identification Task score
#NZwellForm$workerId<- paste0("nz", NZwellForm$workerId)
#length(unique(NZwellForm$word)) #521
# the subset of the participants: select people do not learn Māori at all
#NZwellFormNonMao<-NZwellForm %>% filter(maoriEducation %in% c("0","1"))
#length(unique(NZwellFormNonMao$workerId)) #149
tableEdu<-NZwellFormNonMao%>% dplyr::select(workerId,maoriEducation) %>% distinct_all () %>% group_by(maoriEducation) %>% summarise(n()) # 0: 73, 1:76
# file with 60items to match other data
NZwellForm60<-NZwellFormNonMao[NZwellFormNonMao$word%in%SchoolList$word,]
#length(unique(NZwellForm60$word)) #60
NZwellForm60$dataset<-"adult"
#list(unique(NZwellForm60$word));length(unique(NZwellForm60$workerId)) #149
#names(NZwellForm60)
# original data which contains Identification Task score
OrigWF<-NZwellForm60 %>% dplyr::select (1:29, 60) #%>% dplyr::rename (idt.scr = Exp1.score)
#names(OrigWF);length(unique(OrigWF$workerId)) #149
# school kids data
KidWF<-SchoolWellForm0M%>% dplyr::select (1,3:7, 9, 11:14,16:33,54) %>%
dplyr::rename(nm.nbr = norm.neighbors, p.scr = score.shortv)
#Comparing column names of two files
#fun <- function(OrigWF,KidWF) {
# if(all(names(OrigWF) %in% names(KidWF))) {
# cat('\nNo-Mismatch\n')
# KidWF[names(OrigWF)]
# }
# else {
# cat('\nNo cols : ', setdiff(names(OrigWF), names(KidWF)))
# cat('\nExtra cols : ', setdiff(names(KidWF), names(OrigWF)))
# }
#}
#fun(OrigWF, KidWF) #idt.scr
remove(nb_res)#Finally, we merge all datasets for each task. The resulting datasets are the full datasets used for statistical analysis.
# merge dataset1 and wellForm.
wellFormTask<-rbind(KidWF,OrigWF)
#length(unique(wellFormTask$workerId)) #195
# who is not in the wellform task
#who<-identTask[!identTask$workerId%in%wellFormTask$workerId,]
# remove the participant
#identTask<-identTask[identTask$workerId%in%wellFormTask$workerId,]
remove(KidWF, OrigWF,NZwellForm,SchoolWellForm, SchoolList)
#check numbers of stimuli
#table(identTask$dataset)
#table(wellFormTask$dataset)
#length((unique(identTask$workerId))) #195
#length((unique(identTask$word))) #120
#length((unique(wellFormTask$workerId))) #195
#length((unique(wellFormTask$word))) #60
# stimuli#tableStimuliI<-identTask %>% dplyr::select(word, freq.c, type) %>% distinct_all () %>% group_by(freq.c,type) %>% summarise(n()) # 60 pairs across freq.c
#tableStimuliW<-wellFormTask %>% dplyr::select(word, type,tert_count) %>% distinct_all () %>% group_by(type,tert_count) %>% summarise(n()) # 60 words
identTask<-identTask[identTask$workerId%in%wellFormTask$workerId,]
#length((unique(identTask$workerId))) #1952.4 Overview of Participants’ Demographics
Figure 2 summarizes the distribution of participants. Region was grouped into three categories: south (the South Island), urban (urban areas in the North Island (i.e., Auckland and Wellington)), north (all other areas in the North Island including higher Māori population areas). This categorization is based on participants’ responses to the question “In which region have you spent the largest amount of time since you were 7?”. The majority of our adult participants are non-Māori (92%: European New Zealanders known as Pākehā in New Zealand) and English monolinguals (86%).
#ethnicity
#NZwellFormNonMao %>% dplyr::select(workerId,ethnicity) %>% distinct_all () %>% group_by(ethnicity) %>% summarise(n()) #M:12 (8%), Non-M 136, NS1
#monolinguals
#NZwellFormNonMao %>% dplyr::select(workerId,ethnicity, otherLangs) %>% distinct_all () %>% group_by(ethnicity, otherLangs) %>% summarise(n()) #128 (including 12M) monolinguals.
# extract workerId, age, gender, ethnicity,
origIdentSurvey<- NZwellFormNonMao %>% dplyr::select (1,34,35,36,38,40,48,57) %>% distinct() %>% filter(workerId %in% identTask$workerId) %>% add_column(dataset = "adult") %>% dplyr::rename (age_cat = age) %>% mutate(across('age_cat', str_replace, '60', '60+')) #%>% mutate(age = age %>% as.numeric())
# School
SchSurvey<- SchoolIdent %>% dplyr::select (1,34,35,36,38) %>% distinct() %>% filter(workerId %in% identTask$workerId) %>% mutate(dataset="adolescent")%>% mutate(age_cat = "12-18")
data<-dplyr::bind_rows(origIdentSurvey, SchSurvey)
data<-data %>% mutate(island_cat = case_when(island == "north" ~ "north",
island == "south" ~ "south",
island == "equal" & regionMost %in% c("Canterbury","Nelson Bays","Timaru - Oamaru") ~"south",island == "equal"& regionMost == "Auckland" ~"north") %>% as.factor())
data<-data %>% mutate(region_cat = regionMost,
region_cat = case_when(island_cat == "south" ~ "south",
island == "north" & regionMost %in% c("Otago") ~"south",
regionMost %in% c("Auckland", "Wellington") ~ "urban",
regionMost %in% c("Bay of Plenty","Hawke's Bay","Manawatu","Northland","Taranaki","Wairarapa","Waikato") ~"north") %>% as.factor()) %>% mutate(region_cat=fct_relevel(region_cat,c("south", "urban","north")))
#remove(origIdentSurvey,SchSurvey)
# Population and Dataset Bar Graph
#table(data$dataset)
fig.all<-ggplot(data, aes(x=factor(dataset), fill=factor(dataset)))+
geom_bar(aes(),size=.8,position="dodge",show.legend=FALSE) +
scale_fill_manual(name ="Dataset", values = c("adolescent" = '#619CFF', "adult" = '#F8766D'), guide = guide_legend(reverse = TRUE))+
geom_text(stat='count', aes(label=..count..), vjust=.3, size=3) +
ggtitle("Population")+
ylab('No. of People')
#fig.all
# Age and Dataset Bar Graph
#table(data$age_cat,data$dataset)
fig.age<-ggplot(data, aes(x=factor(age_cat), fill=factor(dataset)))+
geom_bar(aes(),size=.8,position="dodge",show.legend=F) +
scale_x_discrete("Age Group",
labels=c("1"="12-18","2"="18-29","3"="30-39",
"4"="40-49","5"="50-59","6"="60+"))+
scale_fill_manual(name ="dataset", values = c("adolescent" = '#619CFF', "adult" = '#F8766D'), guide = guide_legend(reverse = TRUE))+
geom_text(stat='count', aes(label=..count..), vjust=.3, size=3) +
ggtitle("")+
ylab('No. of People')
#fig.age
#ggplot_build(age)$data
#Gender and Dataset
#table(data$gender, data$dataset)
fig.gender<-ggplot(data, aes(x=factor(gender), fill=factor(dataset)))+
geom_bar(aes(),size=.8,position="dodge",show.legend=F) +
scale_x_discrete("Gender Group",
labels=c("1"="F","2"="M","3"="NB","4"="NS"))+
scale_fill_manual(name ="dataset", values = c("adolescent" = '#619CFF', "adult" = '#F8766D'), guide = guide_legend(reverse = TRUE))+
geom_text(stat='count', aes(label=..count..), position = position_dodge(0.9), size=3) +
ggtitle("")+
ylab('No. of People')
#fig.gender
# Place and Dataset
#table(data$region_cat,data$dataset)
fig.place_region_most <- ggplot(data,aes(x=factor(region_cat), fill=factor(dataset)))+
geom_bar(aes(),size=.8,position="dodge",show.legend=FALSE) +
scale_x_discrete("Region", labels=as.character(c("1"="south","2"="urban","3"="north")))+
scale_fill_manual(name ="dataset", values = c("adolescent" = '#619CFF', "adult" = '#F8766D'), guide = guide_legend(reverse = TRUE))+
geom_text(stat='count', aes(label=..count..), vjust=.3, size=3) +
ggtitle("")+
ylab('No. of People')
#fig.place_within_NZ
ggarrange(fig.all, fig.age, fig.gender, fig.place_region_most)
remove(fig.age, fig.gender, fig.place_region_most, fig.all)
remove(origIdentSurvey,SchSurvey, NZident120, NZwellForm60, NZident0Mao,tableEdu, data)
Figure 2: Overview of Participants’ Demographics
3 Analysis of Adolescent Data (RQ1)
We start by analyzing the adolescent data, in order to address RQ1.
RQ 1. Does exposure to te reo Māori in childhood lead to a proto-lexicon? (i.e., Do school-aged adolescents have a proto-lexicon?)
If the adolescent have a robust proto-lexicon, their ratings in the WRT will correlate significantly with phonotactic score, and in the IDT, they will give significantly higher ratings to words than non-words.
3.1 Age Group
As mentioned above, Adolescent Data contains adolescents whose Māori word knowledge is either 0 or 1 (for a total Māori word knowledge of 1-10). This group is divided into two groups by age (under 15 vs. over 15) in order to examine whether older adolescents showed signs of increased implicit knowledge. Each group has 23 students.
#length(unique(SchoolIdent0M$workerId)) #46
#length(unique(SchoolWellForm0M$workerId)) #46
#a <- SchoolIdent0M[SchoolIdent0M$workerId %in% SchoolWellForm0M$workerId,] # okay
# Age
SchoolIdent0M%>% dplyr::select(workerId,age) %>% distinct() %>% ggplot(aes(x=factor(age), fill=factor(age), label = workerId))+
geom_bar(aes(),size=.8,position="dodge",show.legend=F) +
geom_text(stat='count', aes(label=..count..), vjust=.5) +
ggtitle("Age Group")+ylab('No. of People')
#scale_x_discrete("Age Groups", labels=c("1"="young","2"="old"))+
#ggtitle("Binary Age Group")+ylab('No. of People')
SchoolIdent0M$ageGroup<-"under15"
SchoolIdent0M$ageGroup[which(SchoolIdent0M$age >= "15", )]<-"over15"
#table(SchoolIdent0M$ageGroup) #23 each per group
SchoolWellForm0M$ageGroup<-"under15"
SchoolWellForm0M$ageGroup[which(SchoolWellForm0M$age >= "15", )]<-"over15"
#table(SchoolWellForm0M$ageGroup) #23 each per group
#SchoolWellForm0M %>% dplyr::select(workerId,maoriProf) %>% distinct_all () %>% group_by(maoriProf) %>% summarise(n()) #8:35:3
Figure 3: Age of Students
3.2 Statistical Analyses
The statistical analyses we perform all utilize (logit) ordinal regression, as implemented in the functions clm (fixed-effects only) and clmm (mixed-effects) in the R package ordinal (Christensen 2019). Ordinal regression is used to predict an ordinal dependent variable given one or more independent variables. Likert scales are ordinal scales (i.e., continuous variables) and the scale is not equidistant. Ordinal regression models account for differences in the magnitudes within the thresholds of the scales and predict participant ratings in relation to their predictors, which can be interpreted clearly.
For both tasks, the dependent variable was ‘rating’ on a Likert scale. In the Wellformedness rating task, the scale is a wellformedness judgment rating. In the Word Identification task, the scale is a wordhood confidence rating.
Predictors included linguistic and social/exposure factors, described below in the analysis of each model. Categorical predictors were treatment-, sum- or Helmert-coded as appropriate, and numerical predictors were centered. All models included random intercepts for participant and stimulus, with by-participant random slopes for linguistic predictors and by-stimulus random slopes for social/exposure predictors.
In the fixed-effects setting, we remove interactions one at a time if they do not make a significant contribution to the model’s explanatory power. We then remove main effects if they are not included in any remaining interactions and if they are not statistically significant.
Once there are no more candidate terms for removal, we re-fit the model in a mixed-effects setting, adding random intercepts by participant and item and random slopes by participant for each remaining term. We remove interaction terms from random slopes as required to circumvent problems in running the models. We identify candidate terms for removal if their coefficient is not significant in the model; after performing each removal, we confirm that it was justified via model comparison, using a log-likelihood ratio test.
3.3 Wellformedness Rating Task Ordinal Regression Model
The analysis of the Wellformedness rating task focuses on
participants’ phonotactic knowledge. Models test whether
trends of students’ sensitivity to Māori phonotactics change by age.
The experiment dataset is structured as follows:
Dependent Factor
- enteredResponse is the wellformedness rating for each stimulus.
Independent Factors
- ageGroup (sum-coded) is the classification of each student (over 15 [1] or under 15 [-1]).
- phonotactic score (p.scr) is the phonotactic score for each word. A higher value represents a highly typical Māori (non-)word, while a lower value represents a highly atypical Māori (non-)word.
- neighborhood.occupancy.rate (nm.nbr) is a measure of the number of variations on a stimulus with an edit distance of one that corresponds to a real Māori word and normalized by the maximum number of phonotactically legal neighbors it could potentially have.
Random Effects
- workerId is the unique ID for each participant.
- word is the stimulus used for the rating.
We started a model: c.(p.scr)*ageGroup + c.(nm.nbr).
nonMaoKnowledgeW<- SchoolWellForm0M
#length(unique(nonMaoKnowledgeW$workerId)) #46
#names(nonMaoKnowledgeW)
#str(nonMaoKnowledgeW)
nonMaoKnowledgeW$workerId = as.factor(nonMaoKnowledgeW$workerId)
nonMaoKnowledgeW$enteredResponse = as.factor(nonMaoKnowledgeW$enteredResponse)
nonMaoKnowledgeW$word = as.factor(nonMaoKnowledgeW$word)
nonMaoKnowledgeW$ageGroup = as.factor(nonMaoKnowledgeW$ageGroup)
colnames(nonMaoKnowledgeW)[28] <- "p.scr"
colnames(nonMaoKnowledgeW)[33] <- "nm.nbr"
nonMaoKnowledgeW$ageGroupSum <- as.factor(nonMaoKnowledgeW$ageGroup)
sumC<- contr.sum(2)
colnames(sumC) <- c("_over15") # 1 (over15), -1(under15)
contrasts(nonMaoKnowledgeW$ageGroupSum) = sumC
#contrasts(nonMaoKnowledgeW$ageGroupSum)
# manually coding
#sumCode = matrix(c( 1/2, -1/2), ncol = 1)
#colnames(sumCode) <- c("_over15") # 1 (over15), -1(under15)
#contrasts(nonMaoKnowledgeW$ageGroupSum) = sumCode
#contrasts(nonMaoKnowledgeW$ageGroupSum)
#CHN0W1<- clm(enteredResponse ~ c.(p.scr)*ageGroupSum + c.(nm.nbr), data=nonMaoKnowledgeW)
#summary(CHN0W1) # Two-way interaction is significant - keep all effects.Switch to clmm
#CHN0W2<- clmm(enteredResponse ~ c.(p.scr)*ageGroupSum + c.(nm.nbr) + (1+ c.(p.scr) + c.(nm.nbr)| workerId) + (1+ageGroupSum|word), data=nonMaoKnowledgeW)
#summary(CHN0W2) # interaction is not significant but still phonotactic probability is significant.
#write_rds(CHN0W2, "CHN0W2.rds")
CHN0W2 <- read_rds("docs/CHN0W2.rds")
#count stimuli
#nonMaoKnowledgeW$tert_count <- factor(nonMaoKnowledgeW$tert_count, levels=c("high","medium","low"))
#table(nonMaoKnowledgeW$enteredResponse,nonMaoKnowledgeW$tert_count)
# Get mean ratings for each stimulus per bin
#allBinsScored <- nonMaoKnowledgeW %>%
# group_by(word, bin=tert_count, p.scr) %>% summarise(meanRating =
# mean(as.numeric(enteredResponse))) %>%
# group_by(bin) %>%
# dplyr::summarize(grand_mean = mean(meanRating))The final model is shown in Table 1. The model shows a significant effect of phonotactic score; as seen in the Wellformedness Rating Task Descriptive Statistics, stimuli are rated higher if they have higher phonotactic scores. We also detected a marginal two-way interaction between phonotactic score and ageGroup. This interaction is shown in Figure 4.
names(CHN0W2$coefficients) <- names(CHN0W2$coefficients) %>% str_replace_all("p.scr", "phonotactic score") %>% str_replace_all("nm.nbr", "neighborhood.occupancy.rate")
clm_table(CHN0W2, caption = "Table 1: Ordinal mixed-effects model summary for wellformedness ratings.")| Parameter | Estimate | Std. Error | \(z\) | \(p\) | ||
|---|---|---|---|---|---|---|
| Effects | phonotactic score (centered) | 3.990 | 0.800 | 4.987 | <0.001 | *** |
| ageGroupSum = _over15 | -0.002 | 0.101 | -0.023 | 0.982 | ||
| neighborhood.occupancy.rate (centered) | -2.740 | 1.788 | -1.532 | 0.125 | ||
| phonotactic score (centered) × ageGroupSum = _over15 | 0.968 | 0.569 | 1.701 | 0.089 | . | |
| Thresholds | 1|2 | -2.529 | 0.142 | |||
| 2|3 | -0.710 | 0.132 | ||||
| 3|4 | 0.563 | 0.131 | ||||
| 4|5 | 2.330 | 0.140 |
CHN0W2<- read_rds("docs/CHN0W2.rds")
#post-hoc test
em<-emtrends(CHN0W2, ~ ageGroupSum, var="p.scr",infer =TRUE)
em_table(em) %>%
kable(caption="Table 2: Post hoc Test - Estimated marginal means of linear trends by age group.", label=NA,digits = 3,col.names=c("Frequency","Estimate","Std. Error","lower CL","upper CL","$z$", "$p$", ""),align="lrrrrrrl") %>% kable_styling() | Frequency | Estimate | Std. Error | lower CL | upper CL | \(z\) | \(p\) | |
|---|---|---|---|---|---|---|---|
| over15 | 4.958 | 1.029 | 2.942 | 6.974 | 4.821 | <0.001 | *** |
| under15 | 3.022 | 0.933 | 1.194 | 4.851 | 3.240 | 0.001 | ** |
In order to confirm that the selected model shows no intercorrelation
between predictors, we conducted the Variance Inflation Factor (VIF)
test for each fixed effect. We use a vif() function from the
car package and a polr () function from the
MASS package. First, we fit the ordinal regression model
using the polr (), then pass the fitted model to the vif (). As shown in
Table 3, the VIF test confirms that all fixed effects have a VIF lower
than 5, which indicates that all fixed effects are within an acceptable
range.
CHN0W2 <- read_rds("docs/CHN0W2.rds")
p_CHN0W2<- polr(CHN0W2)
CHN0W2_vif<-vif(p_CHN0W2)
CHN0W2_vif %>%
data.frame() %>%
rownames_to_column(var = "Factor") %>%
setNames(c("Factor", "VIF")) %>%
mutate(
VIF = round(VIF, 2)
) %>%
mutate(Factor = Factor %>%
str_replace_all(.,
"(?<=^|\\:)(?:c\\.\\()([^:]+)(?:\\))(?=$|\\:)",
"\\1 (centered)") %>%
str_replace_all("p.scr", "phonotactic score")%>%
str_replace_all("nm.nbr", "neighborhood.occupancy.rate") %>%
str_replace_all(., fixed(":"), " x "),)%>%
kable(caption="Table 3: Variance Inflation Factor (VIF) for the Wellformedness Rating Task model.", label=NA) %>% kable_styling()| Factor | VIF |
|---|---|
| phonotactic score (centered) | 1.1 |
| ageGroupSum | 1.0 |
| neighborhood.occupancy.rate (centered) | 1.1 |
| phonotactic score (centered) x ageGroupSum | 1.0 |
Figure 4 shows the marginal interaction between ageGroup and phonotactic score. Participants aged over 15 might have even stronger phonotactic intuitions than participants aged under 15.
fig4 <- clm_plotdat(CHN0W2, c("p.scr","ageGroupSum"), xlevels=list(p.scr=25), type="mean") %>%
ggplot(., aes(x=p.scr, y=pred, color=ageGroupSum, fill=ageGroupSum)) +
geom_ribbon(aes(ymin=lci, ymax=uci), alpha=0.3, color=NA) +
geom_line(size=1) +
xlab("Phonotactic Score") +
ylab("Predicted Mean Rating") +
ylim(1, 5) +
scale_x_continuous(guide=guide_axis(check.overlap = TRUE)) +
theme_bw() +
theme(
legend.position="right",
panel.background = element_blank()
)
fig4 + scale_color_manual(name="age group",
values = c("over15" = "blue", "under15" = "red"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="age group",
values = c("over15" = "#619CFF", "under15" = "#F8766D"),
guide = guide_legend(reverse = TRUE))
#remove(modN0W2)
#ggsave("Figure_1A_WRT_Interaction_plot_adolescent.tiff", plot=last_plot(), device = "tiff", scale = 1, dpi = 300)
remove(fig4, CHN0W2)
Figure 4 : Relationship Between Predicted Wellformedness Rating, Phonotactic score and Age Group in the Wellformedness Rating Task Model.
#names(SchoolIdent0M)
#str(SchoolIdent0M)
#change colnames
colnames(SchoolIdent0M)[13] <- "freq.c"
colnames(SchoolIdent0M)[16] <- "n.phon"
colnames(SchoolIdent0M)[28] <- "p.scr"
#names(SchoolIdent0M)
nonMaoKnowledge<-SchoolIdent0M3.4 Identification Task Ordinal Regression Model
The Identification Task (IDT) asks participants to distinguish real words from phonotactically-matched non-words, assessing their proto-lexical knowledge. In order to investigate thoroughly the emergence of Māori proto-lexicon, the adolescent data was analyzed using a (logit) ordinal regression model through the same procedure described in Section 3.2.
The experiment dataset is structured as follows:
Dependent Factor
- enteredResponse is the wordhood confidence rating for each stimulus.
Independent Factors
- type is the type of stimulus (word vs. nonword).
- phonotactic score (p.scr) is the phonotactic score for each word.
- frequency category is the frequency of the word stimuli (mid vs. high).
- ageGroup (sum-coded) is the binary age group (over 15 [1], under 15 [-1]).
Random Effects
- workerId is the unique ID for each participant.
- word is the stimulus used for the rating.
We start with the largest model with a four-way interaction between the independent factors.
# We start with the largest model with a four-way interaction
#names(nonMaoKnowledge)
#length(unique(nonMaoKnowledge$workerId)) #46
nonMaoKnowledge$ageGroup <- as.factor(nonMaoKnowledge$ageGroup)
nonMaoKnowledge$enteredResponse <- as.factor(nonMaoKnowledge$enteredResponse)
nonMaoKnowledge$workerId <- as.factor(nonMaoKnowledge$workerId)
nonMaoKnowledge$freq.c <- factor(nonMaoKnowledge$freq.c ,levels=c("mid","high"))
nonMaoKnowledge$type<- as.factor(nonMaoKnowledge$type)
nonMaoKnowledge$word <- as.factor(nonMaoKnowledge$word)
nonMaoKnowledge$ageGroup <- as.factor(nonMaoKnowledge$ageGroup)
#str(nonMaoKnowledge)
nonMaoKnowledge$ageGroupSum <- as.factor(nonMaoKnowledge$ageGroup)
sumC<- contr.sum(2)
colnames(sumC) <- c("_over15") # 1 (over15), -1(under15)
contrasts(nonMaoKnowledge$ageGroupSum) = sumC
#contrasts(nonMaoKnowledgeW$ageGroupSum)
# 4-way interaction
# We start with a 4-way interaction
#CHN0I1<- clm(enteredResponse ~ c.(p.scr) * type * freq.c * ageGroupSum, data=nonMaoKnowledge)
#summary(CHN0I1) # four-way interaction is not significant - keep all effects. Remove 4-way interaction
#CHN0I2 <- update(CHN0I1, . ~ . -c.(p.scr):type:freq.c:ageGroupSum)
#summary(CHN0I2) # removal candidate :type:freq.c:ageGroupSum
#CHN0I3<- update(CHN0I2, . ~ . -type:freq.c:ageGroupSum)
#summary(CHN0I3) # removal candidate :c.(p.scr):freq.c:ageGroupSum
#CHN0I4 <- update(CHN0I3, . ~ . -c.(p.scr):freq.c:ageGroupSum)
#summary(CHN0I4) # removal candidate :c.(p.scr):type:ageGroupSum
#CHN0I5<- update(CHN0I4, . ~ . -c.(p.scr):type:ageGroupSum)
#summary(CHN0I5) # removal candidate: -type:ageGroupSum
#CHN0I6<- update(CHN0I5, . ~ . -type:ageGroupSum)
#summary(CHN0I6) # removal candidate: freq.c:ageGroupSum
#CHN0I7<- update(CHN0I6, . ~ . -freq.c:ageGroupSum)
#summary(CHN0I7) # no removal candidate, switch to a mixed effects model
# CLMM mixed-effect model
#CHN0I8<-clmm(enteredResponse ~ (c.(p.scr)* type* freq.c) + c.(p.scr)*ageGroupSum+ (1+ageGroupSum|word) + (1+c.(p.scr) + type + freq.c|workerId), data=nonMaoKnowledge)
#summary(CHN0I8) # removal candidate: c.(p.scr):type:freq.c
#CHN0I9<- update(CHN0I8, . ~ . -c.(p.scr):type:freq.c)
#summary(CHN0I9)
#anova(CHN0I8, CHN0I9) # not significant Pr(>Chisq) = 0.429
#CHN0I10<- update(CHN0I9, . ~ . -c.(p.scr):ageGroupSum)
#summary(CHN0I10) # not significant Pr(>Chisq) = 0.2091
#anova(CHN0I9, CHN0I10)
#CHN0I11<- update(CHN0I10, . ~ . -c.(p.scr):type)
#summary(CHN0I11)
#saveRDS(CHN0I11, file ="CHN0I11.rds")
#CHN0I12<-clmm(enteredResponse ~ freq.c * (c.(p.scr) + type) + (1|word) + (1+c.(p.scr) + type + freq.c|workerId), data=nonMaoKnowledge)
#summary(CHN0I12)
#saveRDS(CHN0I12, file = "CHN0I12.rds")
#anova(CHN0I11, CHN0I12)# not significant Pr(>Chisq)= 0.2483
#CHN0I13<-update(CHN0I12, . ~ . -freq.c:type)
#summary(CHN0I13)
#saveRDS(CHN0I13, file = "CHN0I13.rds")
#anova(CHN0I12,CHN0I13) # Pr(>Chisq) = 0.08448, still not significant, take a smaller model
#CHN0I13 <- read_rds("docs/CHN0I13.rds")Table 4 summarizes the best-fitting model. The results of the model show some key findings. Participants were effective in distinguishing words from non-words, as indicated by the significant effect of word type ( b = 0.505 , z = 3.153, p < 0.01), shown in Figure 5. This shows significant evidence of proto-lexical knowledge for non-adult participants. Thus, regarding RQ1, we have shown that adolescents do have protolexical knowledge. In addition, we find a significant interaction between phonotactic score of stimuli and word frequency ( b = −4.925, z = −3.039, p < 0.01), with the effect of phonotactic score being more pronounced in the mid frequency stimuli. However, we did not observe a significant effect of age group in this model.
CHN0I13<- read_rds("docs/CHN0I13.rds")
names(CHN0I13$coefficients) <- names(CHN0I13$coefficients) %>% str_replace_all("p.scr", "phonotactic score")
clm_table(CHN0I13, caption="Table 4: Ordinal mixed-effects model of wordhood confidence ratings.")
#remove(CHN0I13)| Parameter | Estimate | Std. Error | \(z\) | \(p\) | ||
|---|---|---|---|---|---|---|
| Effects | freq.c = high | 0.203 | 0.152 | 1.333 | 0.182 | |
| phonotactic score (centered) | 7.943 | 1.306 | 6.081 | <0.001 | *** | |
| type = word | 0.505 | 0.160 | 3.153 | 0.002 | ** | |
| freq.c = high × phonotactic score (centered) | -4.925 | 1.621 | -3.039 | 0.002 | ** | |
| Thresholds | 1|2 | -2.214 | 0.166 | |||
| 2|3 | -0.537 | 0.162 | ||||
| 3|4 | 0.808 | 0.162 | ||||
| 4|5 | 2.373 | 0.165 |
A post-hoc EMM test (Table 5) confirmed that, despite this interaction, the effect of phonotactic score was still significant even for the high frequency stimulus.
CHN0I13<- read_rds("docs/CHN0I13.rds")
emt<-emtrends(CHN0I13,~freq.c, var="p.scr",infer =TRUE)
em_table(emt) %>%
kable(caption="Table 5: Post hoc test - Estimated marginal means of linear trends by frequency.", label=NA,digits = 3,col.names=c("Frequency","Estimate","Std. Error","lower CL","upper CL","$z$", "$p$", ""),align="lrrrrrrl") %>% kable_styling()
#remove(CHN0I13)| Frequency | Estimate | Std. Error | lower CL | upper CL | \(z\) | \(p\) | |
|---|---|---|---|---|---|---|---|
| mid | 7.943 | 1.306 | 5.383 | 10.502 | 6.081 | <0.001 | *** |
| high | 3.017 | 1.177 | 0.711 | 5.323 | 2.564 | 0.010 | * |
As shown in Table 6, the VIF test confirms that all fixed effects have a VIF lower than 5, which indicates that all fixed effects are within an acceptable range.
CHN0I13<- read_rds("docs/CHN0I13.rds")
p_CHN0I13<- polr(CHN0I13)
CHN0I13_vif<-vif(p_CHN0I13)
CHN0I13_vif %>%
data.frame() %>%
rownames_to_column(var = "Factor") %>%
setNames(c("Factor", "VIF")) %>%
mutate(
VIF = round(VIF, 2)
) %>%
mutate(Factor = Factor %>%
str_replace_all(.,
"(?<=^|\\:)(?:c\\.\\()([^:]+)(?:\\))(?=$|\\:)",
"\\1 (centered)") %>%
str_replace_all("p.scr", "phonotactic score")%>%
str_replace_all("freq.c", "frequency.category") %>%
str_replace_all(., fixed(":"), " x "),)%>%
kable(caption="Table 6: Variance Inflation Factor (VIF) for the Wellformedness Rating Task model.", label=NA) %>% kable_styling()| Factor | VIF |
|---|---|
| frequency.category | 1.00 |
| phonotactic score (centered) | 2.24 |
| type | 1.00 |
| frequency.category x phonotactic score (centered) | 2.24 |
CHN0I13_data <- clm_plotdat(CHN0I13, c("p.scr", "type", "freq.c"), type="mean")
fig5<-CHN0I13_data %>%
ggplot(aes(x = p.scr, y = pred, color = type, fill = type)) +
geom_line(size=1) +
geom_ribbon(aes(ymin = lci, ymax = uci), alpha = 0.3, colour=NA) +
facet_grid(~freq.c) +
#facet_grid(~freq.c,labeller="label_both", switch="y") +
theme_minimal() +
ylim(1, 5) +
theme_bw() +
# theme_minimal() +
xlab("Phonotactic Score") +
ylab("Predicted Mean Rating") +
theme(
panel.grid = element_blank()
)#+facet_grid(. ~ dataset)
fig5 + scale_color_manual(name="type",
values = c("nonword" = "black", "word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="type",
values = c("nonword" = "gray50", "word" = "dodgerblue1"),
guide = guide_legend(reverse = TRUE))
#ggsave("Figure_1B_IDT_Interaction_adolescent.tiff", plot=last_plot(), device = "tiff", scale = 1, dpi = 300)
remove(fig5)
remove(CHN0I13, CHN0I13_data,p_CHN0I13,nonMaoKnowledge,p_CHN0W2, nonMaoKnowledgeW,sumC)
Figure 5: Relationship Between Predicted Wordhood Confidence Rating and Phonotactic score, Stimulus Type, and Frequency category in the Word Identification Task Model.
4 Comparing the Adult and Adolescent Groups (RQ2)
In order to investigate RQ2, we pooled the data from both participant groups, and examined whether there was a significant difference between the adults and the adolescents. In order to compare these groups directly, responses in the adult data that match stimuli in the adolescent data were extracted.
RQ 2. Does continued exposure into adulthood increase the size of the proto-lexicon ? (i.e., Do adults have a larger proto-lexicon than school-aged adolescents?)
4.1 Wellformedness Rating Task Descriptive Statistics
Adolescent Data contains 60 stimuli for the Wellformedness Task whereas Adult Data contains 521 words. We extract responses that match 60 stimuli in Adolescent Data from Adult Data.
Figure 6 shows the mean rating of stimuli in the Wellformedness Rating task by the phonotactic score of the stimuli across two datasets (Adult vs. Adolescent). The results show a positive correlation between the rating of stimuli and their phonotactics score ( r = 0.66 [Adult], r =0.63 [Adolescent]).
#Get mean ratings for each stimulus per bin
allBinsScored <- wellFormTask %>%
group_by(word, bin=tert_count, p.scr,dataset) %>%
dplyr::summarise(
meanRating = mean(as.numeric(enteredResponse))
)
# get boundaries between bins
boundaries <- c(max(allBinsScored$p.scr[allBinsScored$bin == "low"]), max(allBinsScored$p.scr[allBinsScored$bin == "medium"]))
allBinsScored$dataset<-factor(allBinsScored$dataset, levels = c("child","adult"), labels = c("adolescent","adult"))
rating_by_phonotactics <- ggplot(allBinsScored, aes(x=p.scr, y=meanRating, colour = dataset, fill=dataset)) +
geom_point(alpha=0.1, size=4) +
#geom_smooth(method="lm", formula="y~x", alpha=0.7, size=1, color="black") +
geom_smooth(method="lm",aes(x=p.scr, y=meanRating, colour=dataset)) +
geom_vline(xintercept=boundaries, linetype="dotted") +
xlab("Phonotactic score") +
ylab("Mean rating (per stimulus)") +
ylim(1, 5) +
theme_bw() +
theme(
panel.grid = element_blank()
)
rating_by_phonotactics+stat_cor(method="spearman") +scale_color_manual(name="dataset",
values = c("adolescent" = "darkblue", "adult" = "red"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="dataset",
values = c("adolescent" = "#619CFF", "adult" = "#F8766D"),
guide = guide_legend(reverse = TRUE))
remove(allBinsScored, rating_by_phonotactics, boundaries)
Figure 6: Mean Wellformedness Rating by Phonotactic Score. Dotted lines separate items into bins of equal size by phonotactic score.
4.2 Wellformedness Rating Task Ordinal Regression Model
The analysis of the Wellformedness rating task focuses on participants’ phonotactic knowledge. We added neighborhood.occupancy.rate as a control predictor. We predict that non-words with a greater neighborhood density would be more likely to be rated as well formed.
The experiment dataset is structured as follows:
Dependent Factor
- enteredResponse is the wellformedness rating for each stimulus.
Independent Factors
- dataset (sum-coded) is the type of dataset (adult [1], adolescent [-1])
- phonotactic score (p.scr) is the phonotactic score for each word.
- neighborhood.occupancy.rate (nm.nbr) is a measure of the number of variations on a stimulus with an edit distance of one that corresponds to a real Māori word and normalized by the maximum number of phonotactically legal neighbors it could potentially have.
Random Effects
- workerId is the unique ID for each participant.
- word is the stimulus used for the rating.
We started from a model using a 3-way interaction between three independent factors that are described above.
wellFormTask$enteredResponse <- as.factor(wellFormTask$enteredResponse)
wellFormTask$word <- as.factor(wellFormTask$word)
wellFormTask$set<- factor(wellFormTask$dataset)
wellFormTask$workerId<- factor(wellFormTask$workerId)
wellFormTask$dataset.S <- as.factor(wellFormTask$dataset)
sumC<- contr.sum(2)
colnames(sumC) <- c("_adult") # 1 (adult), -1(child) = "adolescent"
contrasts(wellFormTask$dataset.S) = sumC
#contrasts(wellFormTask$dataset.S)
#modW1 <- clm(enteredResponse ~ c.(p.scr) *dataset.S * c.(nm.nbr), data=wellFormTask)
#summary(modW1) # three-way interaction is not significant. Keep all two-way interactions.
#modW2 <- clm(enteredResponse ~ (c.(p.scr) + dataset.S + c.(nm.nbr))^2, data=wellFormTask)
#summary(modW2)# two way interactions significant. add random slopes.
# CLMM
#modW3 <- clmm(enteredResponse ~ (c.(p.scr) + dataset.S + c.(nm.nbr))^2+ (1 + c.(p.scr) + c.(nm.nbr)|workerId) + (1 + dataset.S|word), data=wellFormTask)
#summary(modW3) #
#write_rds(modW3, "modW3.rds")
#modW4 <- update(modW3, . ~ . - dataset.S:c.(nm.nbr))
#summary(modW4) #- c.(p.scr):c.(nm.nbr))
#write_rds(modW4, "modW4.rds")
#anova(modW3,modW4) #Pr(>Chisq) = 0.6662 md4 is better
#modW5<- update(modW4, . ~ . - c.(p.scr):c.(nm.nbr))
#summary(modW5) #
#write_rds(modW5, "modW5.rds")
#anova (modW4, modW5) #Pr(>Chisq) = 0.3177 mod5 is better
#modW6 <- clmm(enteredResponse ~ c.(p.scr) * dataset.S+ (1 + c.(p.scr)|workerId) + (1 +dataset.S|word), data=wellFormTask)
#summary(modW6)
#write_rds(modW6, "modW6.rds")
#anova(modW5,modW6) #Pr(>Chisq) = 0.1281, take a smaller model = mdW6The final model is shown in Table 7. The numeric predictor is centered. The model shows a significant effect of phonotactic score; as seen in the Wellformedness Rating Task Descriptive Statistics, stimuli are rated higher if they have higher phonotactic scores. The results also show a significant interaction between phonotactic score and dataset:adults appear to have stronger phonotactic intuitions than adolescents. This interaction is shown in Figure 7.
modW6 <- readRDS("docs/modW6.rds")
names(modW6$coefficients) <- names(modW6$coefficients) %>% str_replace_all("nm.nbr", "neighbourhood density") %>% str_replace_all("p.scr", "phonotactic score")
clm_table(modW6, caption="Table 7: Ordinal mixed-rffects model summary for wellformedness ratings.", label=NA)| Parameter | Estimate | Std. Error | \(z\) | \(p\) | ||
|---|---|---|---|---|---|---|
| Effects | phonotactic score (centered) | 5.418 | 0.706 | 7.671 | <0.001 | *** |
| dataset.S = _adult | 0.029 | 0.081 | 0.351 | 0.725 | ||
| phonotactic score (centered) × dataset.S = _adult | 1.230 | 0.367 | 3.349 | <0.001 | *** | |
| Thresholds | 1|2 | -2.508 | 0.130 | |||
| 2|3 | -0.712 | 0.124 | ||||
| 3|4 | 0.551 | 0.124 | ||||
| 4|5 | 2.228 | 0.128 |
A post-hoc EMM test (Table 8) confirmed that, despite this interaction, the effect of phonotactic score was still significant even for the adolescent dataset.
#modW6 <- readRDS("docs/modW6.rds")
#emt1<-emtrends(modW6,~dataset.S, var="p.scr",infer =TRUE)
#emt1
tibble(
`Dataset` = c("adult","adolescent"),
`Estimate` = c(6.65, 4.19),
`Std. Error` = c(0.857,0.730),
`lower CL` = c("4.97","2.76"),
`upper CL` = c("8.33"," 5.62"),
`z` = c(7.759, 5.735),
`p.value` = c("<.0001","<.0001"),
) %>%
kbl(caption = "Table 8: Post hoc test - Estimated marginal means of linear trends by dataset.",align="lrrrrrrl") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)| Dataset | Estimate | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|
| adult | 6.65 | 0.857 | 4.97 | 8.33 | 7.759 | <.0001 |
| adolescent | 4.19 | 0.730 | 2.76 | 5.62 | 5.735 | <.0001 |
As shown in Table 9, all fixed effects have a VIF lower than 5, which indicates that all fixed effects are within an acceptable range.
modW6 <- readRDS("docs/modW6.rds")
p_modW6<-polr(modW6)
modW6_vif <-vif(p_modW6)
modW6_vif %>%
data.frame() %>%
rownames_to_column(var = "Factor") %>%
setNames(c("Factor", "VIF")) %>%
mutate(
VIF = round(VIF, 2)
) %>%
mutate(
Factor = Factor %>%
str_replace_all(.,
"(?<=^|\\:)(?:c\\.\\()([^:]+)(?:\\))(?=$|\\:)",
"\\1 (centered)"
) %>%
str_replace_all("p.scr", "phonotactic score")%>%
str_replace_all(., fixed(":"), " x "),
) %>%
kable(caption="Table 9: Variance Inflation Factor (VIF) for the Wellformedness Rating Task model.", label=NA) %>% kable_styling()| Factor | VIF |
|---|---|
| phonotactic score (centered) | 1.05 |
| dataset.S | 1.00 |
| phonotactic score (centered) x dataset.S | 1.05 |
modW6 <- readRDS("docs/modW6.rds")
fig.wellform <- clm_plotdat(modW6, c("p.scr","dataset.S"), xlevels=list(p.scr=25), type="mean") %>%
ggplot(., aes(x=p.scr, y=pred, color=dataset.S, fill=dataset.S)) +
geom_ribbon(aes(ymin=lci, ymax=uci), alpha=0.3, color=NA) +
geom_line(size=1) +
xlab("Phonotactic Score") +
ylab("Predicted Mean Rating") +
ylim(1, 5) +
scale_x_continuous(guide=guide_axis(check.overlap = TRUE)) +
#scale_fill_discrete(name="ponotactic knowledge\n(different sets)") +
theme_bw() +
theme(
legend.position="right",
panel.background = element_blank()
)
#change labels (child to adolescent)
levels(fig.wellform$data$dataset.S)<-c("adult","adolescent")
fig.wellform +scale_color_manual(name="dataset.S",
values = c("adolescent" = "darkblue", "adult" = "red"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="dataset.S",
values = c("adolescent" = "#619CFF", "adult" = "#F8766D"),
guide = guide_legend(reverse = TRUE)) +
theme(text = element_text(size = 14))
#ggsave("Figure_2_WRT_Interaction_with_datasets.tiff", plot=last_plot(), device = "tiff", scale = 1, dpi = 300)
modW6 <- readRDS("docs/modW6.rds")
#em2<-emmip(modW6,dataset.S~p.scr, at=list(dataset.S = c("adult","child"), p.scr= c(-1.2,-1.0,-0.8,-0.6)), CIs=TRUE)
#levels(em2$data$dataset.S)<-c("adult","adolescent")
#levels(em2$data$tvar)<-c("adult","adolescent")
#em2
remove(modW6, p_modW6, modW6_vif, fig.wellform, sumC, wellFormTask,NZwellFormNonMao)
Figure 7: Relationship between Predicted Welformedness Rating, Dataset, and Phonotactic Score.
4.3 Identification Task Descriptive Statistics
The adolescent dataset contains 120 word items for the Identification Task whereas the adult dataset contains 544 word items. We extract responses that match 120 stimuli in the adolescent dataset from the adult datset. As a result, each adult participant responded to 20 to 42 stimuli instead of 120 stimuli.
4.3.1 Mean Wordhood Confidence Ratings
Figure 8 visualizes how the mean confidence rating for a stimulus across datasets differs based on its frequency bin and whether it is a real word or a non-word. Words are generally given higher ratings than non-words, with the size of the difference increasing as the frequency of the word increases. The distinction between ‘Adult Data’ and ‘Adolescent Data’ is apparent in word stimuli. Adults rate words higher than adolescents do but for nonwords stimuli, a rating tendency is less likely to be different.
# Code for Figure S2: Plot mean ratings per bin, per stimulus type
#Get mean ratings for each stimulus per bin
#Get mean ratings for each stimulus per bin
allBins <- list()
bin_names <- c("mid", "high")
for(i in 1:length(bin_names)){
Bin <- identTask[identTask$freq.c==bin_names[i],]
BinMeanRatings <- aggregate(as.numeric(Bin$enteredResponse), by=list(Bin$word, Bin$type,Bin$dataset), mean)
names(BinMeanRatings) <- c("word","type","dataset","meanRating")
BinMeanRatings$bin <- bin_names[i]
allBins[[i]] <- BinMeanRatings
}
remove(i, Bin, BinMeanRatings)
allBins <- do.call(rbind,allBins)
allBins$type = factor(allBins$type)
allBins$dataset = factor(allBins$dataset)
levels(allBins$type) = c("Nonword", "Word")
allBins$bin <- factor(allBins$bin, levels=c("mid","high"))
allBins$dataset<-factor(allBins$dataset, levels = c("child","adult"), labels = c("adolescent","adult"))
# Get grand mean ratings for each stimulus type per bin
binMeans <- allBins %>%
group_by(type, bin, dataset) %>%
dplyr::summarize(grand_mean = mean(meanRating))
ggplot(allBins,aes(x = bin, y = meanRating, fill = type)) +
geom_violin(alpha=0.4, position=position_identity()) +
geom_point(data=binMeans, aes(y=grand_mean, shape=type), size=4) +
geom_line(data=binMeans, aes(x=as.numeric(as.factor(bin)), y=grand_mean, color=type, linetype=type), size=1) +
scale_fill_manual(name="Stimulus type",
values = c("Nonword" = "black", "Word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_shape_manual(name="Stimulus type",
values = c("Nonword" = 21, "Word" = 24),
guide = guide_legend(reverse = TRUE)) +
scale_color_manual(name="Effect estimate",
values = c("Nonword"="black", "Word" = "blue"),
guide="none") +
scale_linetype_manual(name="Effect estimate",
values = c("Nonword" = "solid", "Word" = "dotted"),
guide="none") +
labs(y = "Mean Rating (by stimulus)",x = "Frequency Bin") +
ylim(1, 5) +
theme(axis.text=element_text(size=28),
axis.title=element_text(size=24,face="bold"),
legend.title=element_text(size=22),
legend.text=element_text(size=22)) +
theme_classic()+
facet_grid(. ~ dataset) +
remove(binMeans)
Figure 8: Mean wordhood confidence ratings for each stimulus per frequency category bin across datasets. Points represent mean ratings across all real words and nonwords within each bin.
4.3.2 Mean Rating and Phonotactic Score for each Stimulus per Frequency Bin
Figure 9 expands on the visualization of Figure 8 by adding a dimension for the phonotactic score of the words and non-words within each frequency bin. Across two datasets (top: adolescent, bottom: adult), participants gave higher ratings to words than non-words. Participants also gave higher ratings to both words and non-words with higher phonotactic scores. The effect of phonotactic score was weak in the high-frequency bin (right) compared to the mid-frequency bin (left).
scores <- unique(identTask[,c("word","p.scr")])
allBinsScored <- merge(allBins, scores, by="word")
remove(scores)
discrim_bins_scatter <- ggplot(allBinsScored, aes(x=p.scr, y=meanRating, color=type)) +
geom_point(aes(shape=type), alpha=0.3, size=3) +
geom_smooth(aes(fill=type), method="lm", formula="y~x", alpha=0.6, size=1) +
facet_grid(dataset ~ bin, labeller="label_both", switch="y") +
xlab("Phonotactic score") +
ylab("Mean rating (per stimulus)") +
scale_shape_manual(name="Stimulus type",
values = c("Nonword" = 1, "Word" = 2),
guide = guide_legend(reverse = TRUE)) +
scale_color_manual(name="Stimulus type",
values = c("Nonword"="black", "Word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="Stimulus type",
values = c("Nonword"="gray70", "Word" = "dodgerblue1"),
guide = guide_legend(reverse = TRUE)) +
ylim(1, 5) +
theme_bw() +
theme(
panel.grid = element_blank()
)#+facet_grid(. ~ dataset)
discrim_bins_scatter
# ggsave("plots/discrim_bins_scatter.jpg", plot = discrim_bins_scatter , width = 7, height = 4, dpi = 200)
remove(allBins, allBinsScored, bin_names, discrim_bins_scatter)
Figure 9: Mean wordhood confidence rating vs. phonotactic score for each stimulus for real words and nonwords by frequency category bin across datasets. Lines show correlations within each bin, for each stimulus type.
4.4 Identification Task Ordinal Regression Model
As described in 3.2, in fitting our regression model, we use a
step-down procedure: we start with a large model, and identify
candidates for deletion in a stepwise fashion. For practical reasons, we
first remove as many candidates as possible in a fixed-effects setting
(clm) before moving to a mixed-effects setting
(clmm).
The experiment dataset is structured as follows:
Dependent Factor
- enteredResponse is the wordhood confidence rating for each stimulus.
Independent Factors
- type is the type of stimulus (word vs. nonword).
- phonotactic score (p.scr) is the phonotactic score for each word.
- frequency category is the frequency of the word stimuli (mid vs. high).
- dataset (sum-coded) is the type of dataset (adult [1], adolescent [-1]).
Random Effects
- workerId is the unique ID for each participant.
- word is the stimulus used for the rating.
We ran models starting with a 4-way interaction between four factors:phonotactic score, type, dataset and freq.category
identTask$enteredResponse <- as.factor(identTask$enteredResponse)
identTask$word <- as.factor(identTask$word)
identTask$dataset<- factor(identTask$dataset)
identTask$workerId<- factor(identTask$workerId)
identTask$type<- factor(identTask$type)
#colnames(identTask)[10] <- "freq.c"
identTask$freq.c<- factor(identTask$freq.c)
#names(identTask)
#str(identTask)
#length(unique(identTask$workerId)) #195
identTask$freq.c <- factor(identTask$freq.c , levels=c("mid", "high"))
identTask$dataset.S <- as.factor(identTask$dataset)
sumC<- contr.sum(2)
colnames(sumC) <- c("_adult") # 1 (adult), -1(child) = "adolescent"
contrasts(identTask$dataset.S) = sumC
#contrasts(identTask$dataset.S)
# We start with the largest model with a 4-way interaction without word length
#modI1 <- clm(enteredResponse ~ c.(p.scr) * type *dataset.S *freq.c , data=identTask)
#summary(modI1)
#modI2 <- clm(enteredResponse ~ (c.(p.scr) + type + dataset.S +freq.c)^3 , data=identTask)
#summary(modI2) # three-way interaction is significant - keep all effects. Do a mixed effects model.
# CLMM
#modI3 <- clmm(enteredResponse ~ (c.(p.scr) + type + dataset.S + freq.c)^3 + (1+dataset.S|word) + (1+c.(p.scr) + type + freq.c|workerId), data=identTask)
#summary(modI3)
#write_rds(modI3, "modI3.rds") # no more removal candidate
# three-way interactions are not significant
# removal candidate: c.(p.scr):type:freq.c
#modI4 <- update(modI3, . ~ . -c.(p.scr):type:freq.c)
#summary(modI4)
#write_rds(modI4, "modI4.rds")
# three-way interactions are not significant
# removal candidate: c.(p.scr):type:dataset.S
#anova(modI3,modI4) #Pr(>Chisq) = 0.5188 take a smaller model = modI4
#modI5 <- update(modI4, . ~ . -c.(p.scr):type:dataset.S)
#summary(modI5)
#write_rds(modI5, "modI5.rds")
# three-way interactions are not significant
# removal candidate: type:dataset.S:freq.c
#anova(modI4,modI5) #Pr(>Chisq) = 0.4637 take a smaller model = modI5
#modI6 <- update(modI5, . ~ . -type:dataset.S:freq.c)
#summary(modI6)
#write_rds(modI6, "modI6.rds")
# three-way interactions are not significant
# removal candidate: c.(p.scr):dataset.S:freq.c
#anova(modI5,modI6) #Pr(>Chisq) = 0.2442 take a smaller model = modI6
#modI7 <- update(modI6, . ~ . -c.(p.scr):dataset.S:freq.c)
#summary(modI7)
#write_rds(modI7, "modI7.rds")
# no more three-way interaction
# removal candidate: type:freq.c
#anova(modI6,modI7) #Pr(>Chisq) = 0.1088 take a smaller model = modI7
#modI8 <- update(modI7, . ~ . -type:freq.c)
#summary(modI8)
#write_rds(modI8, "modI8.rds")
# no more three-way interaction
# removal candidate: c.(p.scr):dataset.S
#anova(modI7,modI8) #Pr(>Chisq) = 0.1256 take a smaller model = modI8
#modI9 <- update(modI8, . ~ . -c.(p.scr):dataset.S)
#summary(modI9)
#write_rds(modI9, "modI9.rds")
# no more three-way interaction
# removal candidate: c.(p.scr):type
#anova(modI8,modI9) #Pr(>Chisq) = 0.7813 take a smaller model = modI9
#modI10 <- update(modI9, . ~ . -c.(p.scr):type)
#summary(modI10)
#write_rds(modI10, "modI10.rds")
# no removal candidate
#anova(modI9,modI10) #Pr(>Chisq) = 0.5179 take a smaller model = modI10
#table(identTask$enteredResponse,identTask$type)Table 10 summarizes the best-fitting model. It shows a significant effect of phonotactic score, word-type and word frequency. Two-way interactions are also detected, shown in Figure 11a-c. First, there is a significant interaction between the effect of phonotactic score and frequency category ( b = −4.600, z = −2.761, p < 0.01), with the effect of phonotactic score being more pronounced in mid frequency stimuli. Second, a significant effect of stimuli type (nonword vs. word) interacts with a effect of dataset (adult vs. adolescent) ( b = 0.659, z = 7.216 , p < .001), indicating that adult participants are more able to discriminate real words from non-words than school-aged adolescents.
modI10 <- read_rds("docs/modI10.rds")
names(modI10$coefficients) <- names(modI10$coefficients) %>% str_replace_all("p.scr", "phonotactic score")
clm_table(modI10, caption="Table 10: Ordinal mixed-mffects model of wordhood confidence ratings.")| Parameter | Estimate | Std. Error | \(z\) | \(p\) | ||
|---|---|---|---|---|---|---|
| Effects | phonotactic score (centered) | 7.887 | 1.324 | 5.957 | <0.001 | *** |
| type = word | 1.179 | 0.189 | 6.249 | <0.001 | *** | |
| dataset.S = _adult | -0.103 | 0.094 | -1.098 | 0.272 | ||
| freq.c = high | 0.389 | 0.185 | 2.095 | 0.036 | * | |
| phonotactic score (centered) × freq.c = high | -4.600 | 1.666 | -2.761 | 0.006 | ** | |
| type = word × dataset.S = _adult | 0.659 | 0.091 | 7.216 | <0.001 | *** | |
| dataset.S = _adult × freq.c = high | 0.172 | 0.085 | 2.029 | 0.042 | * | |
| Thresholds | 1|2 | -2.206 | 0.174 | |||
| 2|3 | -0.444 | 0.171 | ||||
| 3|4 | 0.959 | 0.171 | ||||
| 4|5 | 2.536 | 0.173 |
modI10 <- read_rds("docs/modI10.rds")
plot1<-emmip(modI10, type~dataset.S, at=list(type = c("nonword","word"), set= c("adult","child")), CIs=TRUE)
levels(plot1$data$dataset.S)<-c("adult","adolescent")
levels(plot1$data$xvar)<-c("adult","adolescent")
#post=hoc test
#modI10 <- read_rds("docs/modI10.rds")
#emmeans(modI10, revpairwise ~ type|dataset.S,infer=TRUE)
tibble(
`Type` = c("nonword","word"),
`Dataset` = c("adult"),
`Emmean` = c("−0.034", "1.804"),
`Std. Error` = c("0.176","0.181"),
`lower CL` = c("−0.379","1.449"),
`upper CL` = c("0.311","2.159"),
`z` = c("−0.192", "9.969"),
`p.value` = c("0.848", "<.0001"),
) %>%
kbl(caption = "Table 11a-1: Post hoc test - Estimated marginal means for words/nonwords by dataset.",align="llrrrrrr") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)
tibble(
`Type` = c("nonword","word"),
`Dataset` = c("adolescent"),
`Emmean` = c("−0.001", "0.519"),
`Std. Error` = c("0.146","0.152"),
`lower CL` = c("−0.286","0.221"),
`upper CL` = c("0.285","0.817"),
`z` = c("−0.004", "3.417"),
`p.value` = c("0.997", "0.001"),
) %>%
kbl(caption = "",align="llrrrrrr") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)
# contrasts
tibble(
`Contrast` = c("word-nonword"),
`Dataset` = c("adult"),
`Estimate` = c("1.84"),
`Std. Error` = c("0.243"),
`lower CL` = c("1.362"),
`upper CL` = c("2.313"),
`z` = c("7.574"),
`p.value` = c("<.0001"),
) %>%
kbl(caption = "Table 11a-2 Estimated marginal means contrasts for words/nonwords by dataset.",align="llrrrrrr") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)
tibble(
`Contrast` = c("word-nonword"),
`Dataset` = c("adolescent"),
`Estimate` = c(0.52),
`Std. Error` = c(0.170),
`lower CL` = c(0.186),
`upper CL` = c(0.853),
`z` = c(3.052),
`p.value` = c(0.0023),
) %>%
kbl(caption = "",align="llrrrrrr") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)
plot2<-emmip(modI10, freq.c ~ p.scr, at=list(freq.c = c("high","mid"), p.scr= c(-1.2,-1.0,-0.8,-0.6)), CIs=TRUE)
#post-hoc test
#emtrends(modI10,~freq.c, var="p.scr",infer =TRUE)
tibble(
`Frequency` = c("mid","high"),
`Estimate` = c(7.89, 3.29),
`Std. Error` = c(1.32,1.19),
`lower CL` = c(5.292, 0.945),
`upper CL` = c(10.48,5.63),
`z` = c(5.957, 2.751),
`p.value` = c("<.0001","0.0059"),
) %>%
kbl(caption = "Table 11b: Estimated marginal means of linear trends by frequency.",align="lrrrrrrr") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)
plot3<-emmip(modI10, dataset.S~freq.c, at=list(freq.c = c("high","mid"), set= c("adult","child")), CIs=TRUE)
levels(plot3$data$dataset.S)<-c("adult","adolescent")
levels(plot3$data$tvar)<-c("adult","adolescent")
#post-hoc test
#emmeans(modI10, revpairwise ~ freq.c|dataset.S,infer=TRUE)
tibble(
`Frequency` = c("mid","high"),
`Dataset` = c("adult"),
`Emmean` = c(0.605, 1.166),
`Std. Error` = c(0.178,0.177),
`lower CL` = c(0.255, 0.820),
`upper CL` = c(0.954,1.512),
`z` = c(3.388, 6.599),
`p.value` = c("0.0007", "<.0001"),
) %>%
kbl(caption = "Table 11c-1: Estimated marginal means for mid/high by dataset.",align="llrrrrrr") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)
tibble(
`Frequency` = c("mid","high"),
`Dataset` = c("adolescent"),
`Emmean` = c(0.151, 0.367),
`Std. Error` = c(0.145,0.147),
`lower CL` = c("−0.133",0.079),
`upper CL` = c(0.435, 0.655),
`z` = c(1.043, 2.500),
`p.value` = c("0.2971", "0.0124"),
) %>%
kbl(caption = "",align="llrrrrrr") %>%
kable_styling( )%>%
row_spec(1, bold = F, hline_after = T)
tibble(
`Contrast` = c("high-mid"),
`Dataset` = c("adult"),
`Estimate` = c(0.561),
`Std. Error` = c(0.24),
`lower CL` = c(0.091),
`upper CL` = c(1.03),
`z` = c(2.338),
`p.value` = c("0.0194"),
) %>%
kbl(caption = "Table 11c-2: Estimated marginal means contrasts for high/mid by dataset.",align="llrrrrrr") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)
tibble(
`Contrast` = c("high-mid"),
`Dataset` = c("adolescent"),
`Estimate` = c(0.216),
`Std. Error` = c(0.16),
`lower CL` = c("−0.098"),
`upper CL` = c(0.53),
`z` = c(1.350),
`p.value` = c("0.1772"),
) %>%
kbl(caption = "",align="llrrrrrr") %>%
kable_styling()%>%
row_spec(1, bold = F, hline_after = T)| Type | Dataset | Emmean | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|---|
| nonword | adult | −0.034 | 0.176 | −0.379 | 0.311 | −0.192 | 0.848 |
| word | adult | 1.804 | 0.181 | 1.449 | 2.159 | 9.969 | <.0001 |
| Type | Dataset | Emmean | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|---|
| nonword | adolescent | −0.001 | 0.146 | −0.286 | 0.285 | −0.004 | 0.997 |
| word | adolescent | 0.519 | 0.152 | 0.221 | 0.817 | 3.417 | 0.001 |
| Contrast | Dataset | Estimate | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|---|
| word-nonword | adult | 1.84 | 0.243 | 1.362 | 2.313 | 7.574 | <.0001 |
| Contrast | Dataset | Estimate | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|---|
| word-nonword | adolescent | 0.52 | 0.17 | 0.186 | 0.853 | 3.052 | 0.0023 |
| Frequency | Estimate | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|
| mid | 7.89 | 1.32 | 5.292 | 10.48 | 5.957 | <.0001 |
| high | 3.29 | 1.19 | 0.945 | 5.63 | 2.751 | 0.0059 |
| Frequency | Dataset | Emmean | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|---|
| mid | adult | 0.605 | 0.178 | 0.255 | 0.954 | 3.388 | 0.0007 |
| high | adult | 1.166 | 0.177 | 0.820 | 1.512 | 6.599 | <.0001 |
| Frequency | Dataset | Emmean | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|---|
| mid | adolescent | 0.151 | 0.145 | −0.133 | 0.435 | 1.043 | 0.2971 |
| high | adolescent | 0.367 | 0.147 | 0.079 | 0.655 | 2.500 | 0.0124 |
| Contrast | Dataset | Estimate | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|---|
| high-mid | adult | 0.561 | 0.24 | 0.091 | 1.03 | 2.338 | 0.0194 |
| Contrast | Dataset | Estimate | Std. Error | lower CL | upper CL | z | p.value |
|---|---|---|---|---|---|---|---|
| high-mid | adolescent | 0.216 | 0.16 | −0.098 | 0.53 | 1.35 | 0.1772 |
The Variance Inflation Factor (VIF) test in Table 12 confirms that all fixed effects have a VIF lower than 5, which indicates that all fixed effects are within an acceptable range.
modI10 <- read_rds("docs/modI10.rds")
p_modI10<- polr(modI10)
modI10_vif<-vif(p_modI10)
modI10_vif %>%
data.frame() %>%
rownames_to_column(var = "Factor") %>%
setNames(c("Factor", "VIF")) %>%
mutate(
VIF = round(VIF, 2)
) %>%
mutate(
Factor = Factor %>%
str_replace_all(.,
"(?<=^|\\:)(?:c\\.\\()([^:]+)(?:\\))(?=$|\\:)",
"\\1 (centered)"
) %>%
str_replace_all("p.scr", "phonotactic score")%>%
str_replace_all("freq.c", "frequency.category") %>%
str_replace_all(., fixed(":"), " x "),
) %>%
kable(caption="Table 12: Variance Inflation Factor (VIF) for the Identifcaiton Task model.", label=NA) %>% kable_styling()| Factor | VIF |
|---|---|
| phonotactic score (centered) | 2.59 |
| type | 1.02 |
| dataset.S | 3.26 |
| frequency.category | 1.04 |
| phonotactic score (centered) x frequency.category | 2.59 |
| type x dataset.S | 1.90 |
| dataset.S x frequency.category | 2.30 |
Figure 10 shows the Predicted Rating of Stimuli by the three factors shown to interact in the Word Identification Task Model: Phonotactic score, stimulus type and dataset.
modI10_data <- clm_plotdat(modI10, c("p.scr","dataset.S","type"), type="mean") #mutate (type=fct_relevel(type,c("nonword","word")))%>%
fig10 <- modI10_data %>%
ggplot(aes(x = p.scr, y = pred, color = type, fill = type)) +
geom_line(size=1) +
geom_ribbon(aes(ymin = lci, ymax = uci), alpha = 0.3, colour=NA) +
#facet_grid(~dataset.S,labeller="label_both", switch="y") +
facet_grid(~dataset.S) +
theme_minimal() +
ylim(1, 5) +
theme_bw() +
xlab("Phonotactic Score") +
ylab("Predicted Mean Rating") +
theme(
panel.grid = element_blank()
)#+facet_grid(. ~ dataset)
levels(fig10$data$dataset.S)<-c("adult","adolescent")
fig10 + scale_color_manual(name="type",
values = c("nonword" = "black", "word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="type",
values = c("nonword" = "gray50", "word" = "dodgerblue1"),
guide = guide_legend(reverse = TRUE)) +
theme(text = element_text(size = 14))
#ggsave("Figure_3_IDT_Interaction.tiff", plot=last_plot(), device = "tiff", scale = 1, dpi = 300)
Figure 10: Relation Between Predicted Wordhood Confidence Rating, Phonotactic Score, Stimulus Type and Dataset from the Best Model
Figure 11a-c shows the Predicted Rating of Stimuli by the two factors shown to interact significantly in the Word Identification Task Model. We conducted post-hoc EMM tests. (a) stimulus type and dataset: Adults can identify words better than school-aged adolescents (b) frequency category and phonotactic scores: phonotactic scores (= p.scr) is a significant predictor for the confidence rating of mid frequency words (c) dataset and frequency category: The high frequency word is a significant predictor for the confidence rating of adults.
ggarrange(plot1, plot2, plot3)
#ggarrange(plot1, plot2)
remove(plot1, plot2, plot3, modI10, modI10_data, p_modI10, SchoolIdent, SchoolIdent0M,SchoolWellForm0M,sumC,fig10)
Figure 11a-c: Two-Way Interactions in the Word Identification Task Model
5 Effects of Exposure within the Adult Data (RQ3)
In order to investigate RQ3, we used the adult IDT data (n = 149). Because we are no longer comparing directly with the adolescent data, we were able to analyze all available stimuli. This contains 272 word-nonword pairs with three frequency categories (high, mid, low).
RQ 3. If adults possess more knowledge of Māori than adolescents, does degree of exposure to the language play a significant role? (i.e., Do adults with more exposure have a larger protolexicon than adults with less exposure?)
5.1 Adult Data Descriptive Statistics
# original non-Māori-speaking NZers
# IdentTask data
#length(unique(NZidentNonMao$workerId)) #149
NZidentNonMao<-NZidentNonMao[NZidentNonMao$workerId%in%identTask$workerId,] #149
#names(NZidentNonMao)
NZidentNonMao$enteredResponse <- as.factor(NZidentNonMao$enteredResponse)
NZidentNonMao$word <- as.factor(NZidentNonMao$word)
NZidentNonMao$workerId<- factor(NZidentNonMao$workerId)
NZidentNonMao$type<- factor(NZidentNonMao$type)
NZidentNonMao$freq.c<- factor(NZidentNonMao$freq.c)
#colnames(NZidentNonMao)[10] <- "freq.c"
#str(NZidentNonMao)
# stimuli
#tableFreq<-NZidentNonMao%>% dplyr::select(pair.number,freq.c) %>% distinct_all () %>% group_by(freq.c) %>% summarise(n()) # high 64, mid102, low107
# participants
#tableParticipant<-NZidentNonMao%>% dplyr::select(workerId,age,maoriEducation) %>% distinct_all () %>% group_by(maoriEducation,age) %>% summarise(n()) # 0=73, 1=76
# Island Category
#tableIsland<-NZidentNonMao%>% dplyr::select(workerId,island,regionMost) %>% distinct_all () %>% group_by(island) %>% summarise(n()) # equal 5, north 109, south 38
# change age to age_cat and island to island_cat. 5 people who stated "equal" will be classified into "south" or "north" depending on the category "regionMost".
NZidentNonMao<-NZidentNonMao %>% dplyr::rename(age_cat = age) %>% mutate(age_cat = as.factor(age_cat),
island_cat = case_when(island == "north" ~ "north",
island == "south" ~ "south",
island == "equal" & regionMost %in% c("Canterbury","Nelson Bays","Timaru - Oamaru") ~"south",island == "equal"& regionMost == "Auckland" ~"north") %>% as.factor())
#tableIsland<-NZidentNonMao %>% dplyr::select(workerId,island_cat) %>% distinct_all () %>% group_by(island_cat) %>% summarise(n()) #north 108, south 41
#levels(NZidentNonMao$island_cat)
# Frequency Category
# helmert for frequency category
helmertC <- contr.helmert(n = 3)
colnames(helmertC) <- c("2v1", "3v-2") # 1 (low), 2(mid), 3(high)
NZidentNonMao<-NZidentNonMao %>% mutate(freq.c=fct_relevel(freq.c, c("low","mid","high")))
NZidentNonMao$freq.c.HM <- as.factor(NZidentNonMao$freq.c)
contrasts(NZidentNonMao$freq.c.HM) <- helmertC
remove(helmertC)
#contrasts(NZidentNonMao$freq.c.HM)
#levels(NZidentNonMao$freq.c.HM)
#str(NZidentNonMao)
NZidentNonMao <-NZidentNonMao %>% mutate(region_cat = regionMost,
region_cat = case_when(island_cat == "south" ~ "south",
island == "north" & regionMost %in% c("Otago") ~"south",
regionMost %in% c("Auckland", "Wellington") ~ "AKL_WLG",
regionMost %in% c("Bay of Plenty","Hawke's Bay","Manawatu","Northland","Taranaki","Wairarapa","Waikato") ~"north") %>% as.factor()) %>% mutate(region_cat=fct_relevel(region_cat,c("north", "AKL_WLG","south")))Figure 12 shows the mean confidence rating of each stimulus per frequency category and their relation to phonotactic scores across three regional groups in raw data. Although across all frequency category bins, real words receive higher ratings than non-words, there is not much difference between regional categories.
NZidentNonMao<-NZidentNonMao %>% mutate(across('age_cat', str_replace, '60', '60+'))
# Regional categories based on regionMost
#Total 149: Akl 37 Wellington 17 = 54, north = Bay of Plenty 5, Hawke's Bay 3, Manawatu 6, Northland 6, Taranaki 6, Waikato 21,Wairarapa 5 = 52, south = 43
# category based on regionMost
NZidentNonMao<-NZidentNonMao %>% mutate(region_cat = case_when(region_cat == "north" ~ "north",
region_cat == "south" ~ "south",
region_cat == "AKL_WLG" ~ "urban")) %>% mutate(region_cat=fct_relevel(region_cat,c("south", "urban","north")))
tableRegions<-NZidentNonMao %>% dplyr::select(workerId,regionMost,region_cat) %>% distinct_all () %>% group_by(region_cat) %>% summarise(n())
#levels(NZidentNonMao$region_cat)
regionPlot =NZidentNonMao %>%
group_by(region_cat, p.scr, type, freq.c) %>% mutate(freq.c=fct_relevel(freq.c, c("high","mid","low")))%>%
summarise(
mean_response = mean(as.numeric(enteredResponse)),
n = n()
) %>%
ungroup() %>%
ggplot(., aes(x=p.scr, y=mean_response, color=type)) +
geom_point(aes(shape=type), alpha=0.3, size=4) +
geom_smooth(aes(fill=type), method="lm", formula="y~x", alpha=0.6, size=1) +
facet_grid(freq.c ~region_cat, labeller="label_both", switch="y") +
xlab("Phonotactic score") +
ylab("Mean rating (per stimulus)") +
scale_shape_manual(name="Stimulus type ",
#labels = c("pseudo" = "Nonword", "real" = "Word"),
values = c("nonword" = 1, "word" = 2),
guide = guide_legend(reverse = TRUE)) +
scale_color_manual(name="Stimulus type ",
#labels = c("pseudo" = "Nonword", "real" = "Word"),
values = c("nonword" = "black", "word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="Stimulus type ",
#labels = c("pseudo" = "Nonword", "real" = "Word"),
values = c("nonword" = "gray70", "word" = "dodgerblue1"),
guide = guide_legend(reverse = TRUE)) +
theme_bw() +
ggtitle("Effect of Regional Category") +
ylim(1, 5) +
theme(
panel.grid = element_blank(),
)
regionPlot
remove(regionPlot, tableRegions)
Figure 12: Mean confidence rating vs. phonotactic score for each stimulus for real words and nonwords by frequency category bin across there regional groups. Lines show correlations within each bin, for each stimulus type.
Figure 13 (top) shows the mean confidence ratings for each stimulus per frequency category, and their relation to phonotactic scores across five age groups. Figure 14 (bottom) shows the mean confidence rating for each stimulus across age groups on its frequency category. There looked to be some difference between the older group (+50) and younger participants. We ran statistical models on this basis.
agePlot = NZidentNonMao %>%
group_by(age_cat, p.scr, type, freq.c) %>% mutate(freq.c=fct_relevel(freq.c, c("high","mid","low")))%>%
summarise(
mean_response = mean(as.numeric(enteredResponse)),
n = n()
) %>%
ungroup() %>%
ggplot(., aes(x=p.scr, y=mean_response, color=type)) +
geom_point(aes(shape=type), alpha=0.3, size=4) +
geom_smooth(aes(fill=type), method="lm", formula="y~x", alpha=0.6, size=1) +
facet_grid(freq.c ~age_cat, labeller="label_both", switch="y") +
xlab("Phonotactic score") +
ylab("Mean rating (per stimulus)") +
scale_shape_manual(name="Stimulus type ",
#labels = c("pseudo" = "Nonword", "real" = "Word"),
values = c("nonword" = 1, "word" = 2),
guide = guide_legend(reverse = TRUE)) +
scale_color_manual(name="Stimulus type ",
#labels = c("pseudo" = "Nonword", "real" = "Word"),
values = c("nonword" = "black", "word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="Stimulus type ",
#labels = c("pseudo" = "Nonword", "real" = "Word"),
values = c("nonword" = "gray70", "word" = "dodgerblue1"),
guide = guide_legend(reverse = TRUE)) +
theme_bw() +
ggtitle("Effect of Age Category") +
ylim(1, 5) +
theme(
panel.grid = element_blank(),
)
agePlot
#Get mean ratings for each stimulus per bin
allBins <- list()
bin_names <- c("low","mid", "high")
for(i in 1:length(bin_names)){
Bin <- NZidentNonMao[NZidentNonMao$freq.c==bin_names[i],]
BinMeanRatings <- aggregate(as.numeric(Bin$enteredResponse), by=list(Bin$word, Bin$type, Bin$age_cat), mean)
names(BinMeanRatings) <- c("word","type","age_cat","meanRating")
BinMeanRatings$bin <- bin_names[i]
allBins[[i]] <- BinMeanRatings
}
remove(i, Bin, BinMeanRatings)
allBins <- do.call(rbind,allBins)
allBins$type = factor(allBins$type)
allBins$set = factor(allBins$age_cat)
levels(allBins$type) = c("Nonword", "Word")
allBins$bin <- factor(allBins$bin, levels=c("low","mid","high"))
# Get grand mean ratings for each stimulus type per bin
binMeans <- allBins %>%
group_by(type, bin,age_cat) %>%
dplyr::summarize(grand_mean = mean(meanRating))
ggplot(allBins,aes(x = age_cat, y = meanRating, fill = type)) +
geom_violin(alpha=0.4, position=position_identity()) +
geom_point(data=binMeans, aes(y=grand_mean, shape=type), size=4) +
geom_line(data=binMeans, aes(x=as.numeric(as.factor(age_cat)), y=grand_mean, color=type, linetype=type), size=1) +
scale_fill_manual(name="Stimulus type",
values = c("Nonword" = "black", "Word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_shape_manual(name="Stimulus type",
values = c("Nonword" = 21, "Word" = 24),
guide = guide_legend(reverse = TRUE)) +
scale_color_manual(name="Effect estimate",
values = c("Nonword"="black", "Word" = "blue"),
guide="none") +
scale_linetype_manual(name="Effect estimate",
values = c("Nonword" = "solid", "Word" = "dotted"),
guide="none") +
labs(y = "Mean Rating (by stimulus)",x = "Frequency Bin") +
ylim(1, 5) +
theme(axis.text=element_text(size=28),
axis.title=element_text(size=24,face="bold"),
legend.title=element_text(size=22),
legend.text=element_text(size=22)) +
theme_classic()+
facet_grid(. ~ bin) +
remove(binMeans, allBins,allBinsScored, agePlot)
# age group
tableAge<-NZidentNonMao %>% dplyr::select(workerId,age_cat) %>% distinct_all () %>% group_by(age_cat) %>% summarise(n()) # "18-29" 60 : "30-39" 28': "40-49" 22:"50-59" 19 :60" 20
NZidentNonMao<-NZidentNonMao %>% mutate(age_B = age_cat,
age_B = case_when(age_cat %in% c("18-29", "30-39", "40-49") ~ "young",
age_cat %in% c("50-59", "60+") ~"old") %>% as.factor()) %>% mutate(age_B=fct_relevel(age_B, c("young","old")))
#levels(NZidentNonMao$age_B) #110,39
remove(tableAge)
Figure 13 (top) and Figure 14 (bottom): Mean confidence rating vs. phonotactic score for each stimulus for real words and nonwords by frequency category bin. Lines show correlations within each bin, for each stimulus type (top). Mean confidence rating for each stimulus for real words and nonwords per frequency category bin (bottom).
5.2 Overview of Participants’ Demographics
Figure 15 summarizes the distribution of participants.
# Age and Māori Education Bar Graph
fig.education<-NZidentNonMao%>% dplyr::select(workerId,age_cat,maoriEducation) %>% distinct() %>%
ggplot(aes(x=factor(maoriEducation), fill=factor(age_cat)))+
geom_bar(aes(),size=.8,position="dodge",show.legend=T) +
scale_x_discrete("Maori Education",
labels=c("0"="never","1"="at school"))+
geom_text(stat='count', aes(label=..count..), position = position_dodge(0.9), size=3) +
ggtitle("Eduation by Age")+
ylab('No. of People')
#fig.education
# Binary age group
fig.age2<-NZidentNonMao %>% dplyr::select(workerId,age_B) %>% distinct() %>% ggplot(aes(x=factor(age_B), fill=factor(age_B), label = workerId))+
geom_bar(aes(),size=.8,position="dodge",show.legend=F) +
geom_text(stat='count', aes(label=..count..), vjust=.5) +
scale_x_discrete("Age Groups",
labels=c("1"="young","2"="old"))+
ggtitle("Binary Age Group")+ylab('No. of People')
#fig.age2 #110,39
#Regional group
fig.region<-NZidentNonMao %>% dplyr::select(workerId,region_cat) %>% distinct() %>%
ggplot(aes(x=factor(region_cat), fill=factor(region_cat)))+
geom_bar(aes(),size=.8,position="dodge",show.legend=F) +
geom_text(stat='count', aes(label=..count..), vjust=.5) +
scale_x_discrete("Region Group",
labels=c("1"="north","2"="urban","3"="south"))+
ggtitle("Reginal Groups")+ylab('No. of People')
#fig.region
# Exposure
#table(data$island,data$dataset)
#cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#7CAE00")
fig.maori.exposure<-NZidentNonMao %>% dplyr::select(workerId,maoriExpo) %>% distinct() %>%
ggplot(aes(x=factor(maoriExpo), fill=factor(maoriExpo)))+
geom_bar(aes(), size =.1, show.legend = F) +
scale_x_discrete("Exposure to Māori",labels=as.character(c("2"=2,"3"=3,"4"=4,"5"=5, "6"=6,"7"=7,"8"=8,"9"=9,"10"=10)))+
geom_text(stat='count', aes(label=..count..), vjust=1) +
#scale_fill_manual(name ="workerId", values = cbPalette, guide = guide_legend(reverse = TRUE))+
ggtitle("Exposure to Māori")+
ylab('No. of People')
#fig.maori.exposure
ggarrange(fig.age2, fig.region, fig.maori.exposure, fig.education)
remove(fig.age2,fig.education,fig.maori.exposure,fig.region)
#sd(NZidentNonMao$maoriExpo) #2.002554
#summary(NZidentNonMao$maoriExpo)#7.10
Figure 15: Overview of participants demographics
5.3 Identification Task Ordinal Regression Model
We focus on three variables which may be linked to the degree of exposure; age, region, and frequency of exposure to Māori based on the self-reported questionnaire.
The experiment dataset is structured as follows:
Dependent Factor
- enteredResponse is the wordhood confidence rating for each stimulus.
Independent Factors
- type is the type of stimulus (word vs. nonword).
- phonotactic score (p.scr) is the phonotactic score for each word.
- frequency category is the frequency of the word stimuli. Helmert contrasts between low vs.mid, and low & mid vs. high.
- age_B is the binary age group, either over 50 or under 50.
- region_cat is each NZ-based participant’s place where they
have spent the largest amount of time since they were 7. (categorized
into three classification,
south= the South Island,urban= urban areas in the North Island (i.e., Auckland, Wellington),north= other areas in the North Island including higher Māori population ares in the North Island) - exposure is each participant’s level of exposure to Māori (with a log-like scale ranging from 2 to 10, combining exposure via media and in person).
Random Effects
- workerId is the unique ID for each participant.
- word is the stimulus used for the rating.
We started from a model using 3-way interactions between five factors: (enteredResponse ~ (c.(p.scr)* type)(age_B+region_cat+maoriExpo) + (typefreq.c.HM)*(age_B+region_cat+maoriExpo). Participant and stimulus were random effects, with linguistic factors as a random slope for ‘participant’ and sociolinguistic factors as a random slope for ‘stimulus’.
#length(unique(NZidentNonMao$workerId)) =149
#Age (binary), Region (three categories), Māori exposure (numeric)
#str(NZidentNonMao)
#levels(NZidentNonMao$region_cat) #south-urban-north
#INM0_sb<- clm(enteredResponse ~ (c.(p.scr)*type)*(age_B+region_cat+c.(maoriExpo)) + (type*freq.c.HM)*(age_B+region_cat+c.(maoriExpo)), data=NZidentNonMao)
#summary(INM0_sb) # removal candidate: type:region_cat:freq.c.HM
#INM1_sb<- update(INM0_sb, . ~ . - type:region_cat:freq.c.HM)
#summary(INM1_sb) # removal candidate:c.(p.scr):type:age_B
#INM2_sb<- update(INM1_sb, . ~ . - c.(p.scr):type:age_B)
#summary(INM2_sb) # removal candidate: c.(p.scr):typeword:c.(maoriExpo)
#INM3_sb<- update(INM2_sb, . ~ . - c.(p.scr):type:c.(maoriExpo))
#summary(INM3_sb) # removal candidate: c.(p.scr):type:region_cat
#INM4_sb<- update(INM3_sb, . ~ . - c.(p.scr):type:region_cat)
#summary(INM4_sb) # removal candidate:c.(p.scr):age_B
#INM5_sb<- update(INM4_sb, . ~ . - c.(p.scr):age_B)
#summary(INM5_sb) # removal candidate: type:region_cat
#INM6_sb<- update(INM5_sb, . ~ . - type:region_cat)
#summary(INM6_sb) # removal candidate: c.(p.scr):maoriExpo
#INM7_sb<- update(INM6_sb, . ~ . - c.(p.scr):c.(maoriExpo))
#summary(INM7_sb) # removal candidate: type:age_B:freq.c.HM
#INM8_sb<- update(INM7_sb, . ~ . - type:age_B:freq.c.HM)
#summary(INM8_sb) # removal candidate:age_B:freq.c.HM
#INM9_sb<- update(INM8_sb, . ~ . - age_B:freq.c.HM)
#summary(INM9_sb) # removal candidate:region_cat:freq.c.HM
#INM10_sb<- update(INM9_sb, . ~ . - region_cat:freq.c.HM)
#summary(INM10_sb) # removal candidate:Switch to mixed-effects
## CLMM
#INM11_sb<- clmm(enteredResponse ~ type*c.(maoriExpo)*freq.c.HM+c.(p.scr)*(type+region_cat) + type*age_B + (1 + c.(p.scr)+type+freq.c.HM|workerId) + (1 + age_B+region_cat+c.(maoriExpo)|word), data=NZidentNonMao)
#summary(INM11_sb)
#write_rds(INM11_sb,"INM11_sb.rds")
#INM12_sb<- clmm(enteredResponse ~ type*c.(maoriExpo)*freq.c.HM+c.(p.scr)*region_cat + type*age_B + (1 + c.(p.scr)+type+freq.c.HM|workerId) + (1 + age_B+region_cat+c.(maoriExpo)|word), data=NZidentNonMao)
#summary(INM12_sb)
#write_rds(INM12_sb,"INM12_sb.rds")
#anova(INM11_sb, INM12_sb) #Pr(>Chisq) = 0.2319, take a smaler model, INM12_sbTable 13 summarizes the best-fitting model. All numeric predictors are centered, and frequency category is Helmert-coded. We find a three-way interaction between the effects of word, word frequency and exposure level ( b = 0.030 , z = 2.806 , p < 0.01). Participants with higher exposure to Māori rate words higher. A post-hoc EMM test confirms that the real word in the high frequency bin is likely to be rated higher according to participants’ level of exposure to Māori (freq.c = High: word b = 0.1178, z = 2.989, p = 0.0028. A significant two-way interaction between the effects of phonotactic score and region: north is also detected ( b = 1.359, z = 2.330, p = 0.020), indicating participants in areas with many Māori speakers are more sensitive to phonotactics than participants in areas with few Māori speakers. This is in line with our prediction. The effects of age and word-type are also interacted ( b = 0.206, z = 2.301, p = 0.021), indicating participants aged over 50 can identify real words better than younger participants. A post-hoc EMM test confirms that this is a significant interaction (type = word: old - young, b = 0.284, z = 1.987, p = 0.047). These interactions are shown in Figure 16 and 17a-b.
Note that we ran another model in which the level of exposure to Māori was categorized binary. The trend of results is the same as our final model which is shown in Table 13.
INM12_sb<- read_rds("docs/INM12_sb.rds")
names(INM12_sb$coefficients) <- names(INM12_sb$coefficients) %>% str_replace_all("p.scr", "phonotactic score")
clm_table(INM12_sb, caption="Table 13: Ordinal mixed-effects model of wordhood confidence ratings.")
#remove(INM12_sb)| Parameter | Estimate | Std. Error | \(z\) | \(p\) | ||
|---|---|---|---|---|---|---|
| Effects | type = word | 0.887 | 0.105 | 8.434 | <0.001 | *** |
| maoriExpo (centered) | 0.035 | 0.032 | 1.088 | 0.277 | ||
| freq.c.HM = 2v1 | 0.058 | 0.078 | 0.743 | 0.458 | ||
| freq.c.HM = 3v-2 | 0.085 | 0.054 | 1.581 | 0.114 | ||
| phonotactic score (centered) | 3.696 | 0.581 | 6.356 | <0.001 | *** | |
| region_cat = urban | 0.019 | 0.148 | 0.131 | 0.896 | ||
| region_cat = north | 0.132 | 0.150 | 0.881 | 0.378 | ||
| age_B = old | 0.077 | 0.148 | 0.522 | 0.601 | ||
| type = word × maoriExpo (centered) | 0.029 | 0.022 | 1.325 | 0.185 | ||
| type = word × freq.c.HM = 2v1 | 0.119 | 0.111 | 1.077 | 0.281 | ||
| type = word × freq.c.HM = 3v-2 | 0.454 | 0.074 | 6.103 | <0.001 | *** | |
| maoriExpo (centered) × freq.c.HM = 2v1 | -0.005 | 0.013 | -0.407 | 0.684 | ||
| maoriExpo (centered) × freq.c.HM = 3v-2 | -0.003 | 0.009 | -0.375 | 0.708 | ||
| phonotactic score (centered) × region_cat = urban | 0.680 | 0.564 | 1.206 | 0.228 | ||
| phonotactic score (centered) × region_cat = north | 1.359 | 0.583 | 2.330 | 0.020 | * | |
| type = word × age_B = old | 0.206 | 0.090 | 2.301 | 0.021 | * | |
| type = word × maoriExpo (centered) × freq.c.HM = 2v1 | 0.004 | 0.018 | 0.230 | 0.818 | ||
| type = word × maoriExpo (centered) × freq.c.HM = 3v-2 | 0.030 | 0.011 | 2.806 | 0.005 | ** | |
| Thresholds | 1|2 | -2.610 | 0.138 | |||
| 2|3 | -0.635 | 0.136 | ||||
| 3|4 | 0.935 | 0.136 | ||||
| 4|5 | 2.524 | 0.137 |
INM12_sb<- read_rds("docs/INM12_sb.rds")
# posthoc test for the interaction between type, exposure level to Māori, and frequency categories.
#emmeans(INM12_sb,~type*maoriExpo*freq.c.HM, adjust="none") %>% pairs(simple = "type", reverse=TRUE, infer=TRUE)
emt1<-emtrends(INM12_sb,~ type|freq.c.HM, var="maoriExpo",adjust="none", infer =TRUE)
# post=hoc test frequency
# Table 14a
summary(emt1)%>% mutate(across(where(is.numeric), round, 3)) %>% mutate(
significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "\\*",
p.value < 0.1 ~ ".",
TRUE ~ ""),
p.value = ifelse(p.value<0.001, "<0.001", format(round(p.value, 3), nsmall=3))) %>%
dplyr::select(-df) %>%
kable(digits = 3,col.names=c("Type","Frequency","Estimate","Std. Error","lower CL", "upper CL", "$z$", "$p$", ""), align="llrrrrrr",
caption = "Table 14a. Estimated marginal means for words/nonwords by frequency.") %>% kable_styling()
# postdoc test for the interaction between phonotactic score and region categories
plot4<-emmip(INM12_sb,region_cat~p.scr, at=list (p.scr = c(-1.2,-1.0,-0.8,-0.6)), CIs=TRUE)
emt2<-emtrends(INM12_sb,pairwise ~region_cat, var="p.scr", adjust="none", infer = TRUE)
# Table 14b-1
summary(emt2$emtrends)%>% mutate(across(where(is.numeric), round, 3)) %>% mutate(
significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "\\*",
p.value < 0.1 ~ ".",
TRUE ~ ""),
p.value = ifelse(p.value<0.001, "<0.001", format(round(p.value, 3), nsmall=3))) %>%
dplyr::select(-df) %>%
kable(digits = 3,col.names=c("Region","Estimate","Std. Error","lower CL", "upper CL", "$z$", "$p$", ""), align="lrrrrrrr",
caption = "Table 14b-1. Estimated marginal means of linear trends by region.") %>% kable_styling()
# Table 14b-2
summary(emt2$contrasts)%>% mutate(across(where(is.numeric), round, 3)) %>% mutate(
significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "\\*",
p.value < 0.1 ~ ".",
TRUE ~ ""),
p.value = ifelse(p.value<0.001, "<0.001", format(round(p.value, 3), nsmall=3))) %>%
dplyr::select(-df) %>%
kable(digits = 3,col.names=c("Contrast","Estimate","Std. Error","lower CL", "upper CL", "$z$", "$p$", ""), align="lrrrrrrr",
caption = "Table 14b-2. Estimated marginal means contrasts for phonotactic scores by region.") %>% kable_styling()
# postdoc test for the interaction between word types score and age categories
plot5<-emmip(INM12_sb,type~age_B, CIs=TRUE)
emm_t.a<-emmeans(INM12_sb, revpairwise ~ type|age_B)
# Table 14c
summary(emm_t.a$emmeans) %>% dplyr::select(-df) %>%
kable(digits = 3,col.names=c("Type","Age","Emmean","Std. Error","lower CL", "upper CL"), align="llrrrr",
caption = "Table 14c. Estimated marginal means for words/nonwords by age.") %>% kable_styling()
#Table 14e: contrasts
emm_a.t<-emmeans(INM12_sb, revpairwise ~ age_B|type)
summary(emm_a.t$contrasts) %>% mutate(
significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "\\*",
p.value < 0.1 ~ ".",
TRUE ~ ""),
p.value = ifelse(p.value<0.001, "<0.001", format(round(p.value, 3), nsmall=3))) %>%
dplyr::select(-df) %>%
kable(digits = 3,col.names=c("Contrast","Type","Estimate","Std. Error","$z$", "$p$", ""), align="llrrrrl",
caption = "Table 14e. Estimated marginal means contrasts for old/young by type.") %>% kable_styling()| Type | Frequency | Estimate | Std. Error | lower CL | upper CL | \(z\) | \(p\) | |
|---|---|---|---|---|---|---|---|---|
| nonword | low | 0.044 | 0.038 | -0.030 | 0.117 | 1.166 | 0.244 | |
| word | low | 0.038 | 0.033 | -0.027 | 0.103 | 1.138 | 0.255 | |
| nonword | mid | 0.033 | 0.037 | -0.038 | 0.105 | 0.912 | 0.362 | |
| word | mid | 0.036 | 0.034 | -0.032 | 0.103 | 1.042 | 0.297 | |
| nonword | high | 0.028 | 0.035 | -0.040 | 0.096 | 0.814 | 0.416 | |
| word | high | 0.118 | 0.039 | 0.041 | 0.195 | 2.989 | 0.003 | ** |
| Region | Estimate | Std. Error | lower CL | upper CL | \(z\) | \(p\) | |
|---|---|---|---|---|---|---|---|
| south | 3.696 | 0.581 | 2.556 | 4.835 | 6.356 | <0.001 | *** |
| urban | 4.376 | 0.605 | 3.190 | 5.562 | 7.232 | <0.001 | *** |
| north | 5.055 | 0.629 | 3.822 | 6.287 | 8.039 | <0.001 | *** |
| Contrast | Estimate | Std. Error | lower CL | upper CL | \(z\) | \(p\) | |
|---|---|---|---|---|---|---|---|
| south - urban | -0.680 | 0.564 | -1.786 | 0.425 | -1.206 | 0.228 | |
| south - north | -1.359 | 0.583 | -2.502 | -0.216 | -2.330 | 0.020 | * |
| urban - north | -0.679 | 0.546 | -1.749 | 0.392 | -1.242 | 0.214 |
| Type | Age | Emmean | Std. Error | lower CL | upper CL |
|---|---|---|---|---|---|
| nonword | young | -0.003 | 0.100 | -0.199 | 0.193 |
| word | young | 0.884 | 0.099 | 0.690 | 1.078 |
| nonword | old | 0.074 | 0.147 | -0.214 | 0.363 |
| word | old | 1.168 | 0.144 | 0.886 | 1.449 |
| Contrast | Type | Estimate | Std. Error | \(z\) | \(p\) | |
|---|---|---|---|---|---|---|
| old - young | nonword | 0.077 | 0.148 | 0.522 | 0.601 | |
| old - young | word | 0.284 | 0.143 | 1.987 | 0.047 | * |
In Table 15, we report Generalized Variance Inflation Factors (GVIFs)
for the model. Same as other models, we used the function polr () within
the package MASS and the function vif () within the package
car. This produces GVIF ,DF (degrees of freedom), and
‘GVIF^(1/(2*DF))’ (GVIF is normalized by the number of DF of GVIF) for
each factor. We use the square-root of GVIF (Fox and Monette 1992)
The VIF test confirms that all fixed effects have a GVIF (i.e., converted GVIF) lower than 5, which indicates that all fixed effects are within an acceptable range.
INM12_sb<- read_rds("docs/INM12_sb.rds")
p_INM12_sb<-polr(INM12_sb)
INM12_sb_vif<-vif(p_INM12_sb)
INM12_sb_vif%>%
data.frame() %>%
setNames(c("GVIF", "DF", "Converted_GVIF")) %>%
rownames_to_column(var = "Factor") %>%
mutate(GVIF = round(GVIF, 2)) %>%
mutate(
Factor = Factor %>%
str_replace_all(.,
"(?<=^|\\:)(?:c\\.\\()([^:]+)(?:\\))(?=$|\\:)",
"\\1 (centered)"
) %>%
str_replace_all("p.scr", "phonotactic score")%>%
str_replace_all("freq.c.HM", "frequency.category") %>%
str_replace_all(., fixed(":"), " x "),
Converted_GVIF = Converted_GVIF ^ 2
) %>%
mutate(Converted_GVIF= round(Converted_GVIF, 2)) %>%
kable(caption="Table 15: Generalised Variance Inflation Factor (GVIF) for the Word Identification Task model.", label=NA) %>%
kable_styling()| Factor | GVIF | DF | Converted_GVIF |
|---|---|---|---|
| type | 1.35 | 1 | 1.35 |
| maoriExpo (centered) | 1.95 | 1 | 1.95 |
| frequency.category | 3.56 | 2 | 1.89 |
| phonotactic score (centered) | 3.20 | 1 | 3.20 |
| region_cat | 1.04 | 2 | 1.02 |
| age_B | 1.95 | 1 | 1.95 |
| type x maoriExpo (centered) | 1.93 | 1 | 1.93 |
| type x frequency.category | 3.57 | 2 | 1.89 |
| maoriExpo (centered) x frequency.category | 3.66 | 2 | 1.91 |
| phonotactic score (centered) x region_cat | 3.20 | 2 | 1.79 |
| type x age_B | 2.26 | 1 | 2.26 |
| type x maoriExpo (centered) x frequency.category | 3.67 | 2 | 1.91 |
Figure 16 shows the predicted rating of stimuli by the four factors shown to interact significantly in the Identification Task model.
INM12_sb<- read_rds("docs/INM12_sb.rds")
INM12_sb_data <- clm_plotdat(INM12_sb, c("type", "maoriExpo", "freq.c.HM"), type="mean")
fig14<-INM12_sb_data %>%
ggplot(aes(x = maoriExpo, y = pred, color = type, fill = type)) +
geom_line(size = 1) +
geom_ribbon(aes(ymin = lci, ymax = uci), alpha = 0.3, colour=NA) +
facet_grid(~freq.c.HM) +
#theme_minimal() +
theme_bw() +
xlab("Maori Exposure") +
ylab("Predicted Mean Rating") +
theme(
panel.grid = element_blank()
)+
scale_y_continuous(limits=c(0, 5), breaks=c(1, 2, 3, 4, 5))
fig14 + scale_color_manual(name="type",
values = c("nonword" = "black", "word" = "blue"),
guide = guide_legend(reverse = TRUE)) +
scale_fill_manual(name="type",
values = c("nonword" = "gray50", "word" = "dodgerblue1"),
guide = guide_legend(reverse = TRUE)) +
theme(text = element_text(size = 14))
#ggsave("Figure_4_IDT_Interactions_adult_plot.tiff", plot=last_plot(), device = "tiff", scale = 1, dpi = 300)
Figure 16 : Relationship Between Phonotactic Score, Stimulus Type, Frequency Category Bin and Dataset from the Final Model
Figure 17a-b shows the predicted rating of stimuli by the two factors
shown to interact significantly in the Identification Task final model.
Post-hoc EMM tests confirm the following results. (a) phonotactic score
and regional category: participants in higher Māori density area in the
North Island (i.e., north) are more sensitive to Māori
phonotactics than participants in the South Island (i.e.,
south) where there are notably fewer Māori speakers (b)
word type category and age group:participants aged over 50 can identify
real words better than younger participants. Despite this interaction,
the effect of word type was still significant even for the younger
people.
ggarrange(plot4, plot5)
remove(INM12_sb, INM12_sb_data, INM12_sb_vif,p_INM12_sb, plot4, plot5,fig14, identTask, NZidentNonMao)
Figure 17a-b: Two Two-Way Interactions
Reference
Fox, J., & Monette, G. (1992). Generalized collinearity diagnostics. Journal of the American Statistical Association, 87(417), 178-183.
Panther, F. A., Mattingley, W., Todd, S., Hay, J. & King, J.(2023). Proto-Lexicon Size and Phonotactic Knowledge are Linked in Non-Māori Speaking New Zealand Adults. Laboratory Phonology 14(1).