F1 F0 cue-weighting

1. Load necessary materials 📚

1.1. Load necessary libraries

library(metR)        # Tools for meteorological data processing

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
There were 50 or more warnings (use warnings() to see the first 50)

library(readr)       # For reading CSV files efficiently
library(dplyr)       # For data manipulation and transformation

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(tidyr)       # For tidying data (reshaping)
library(ggplot2)     # For creating complex data visualizations
library(lme4)        # For fitting linear mixed-effects models

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

library(emmeans)     # For computing estimated marginal means (post hoc analysis)

Warning: package ‘emmeans’ was built under R version 4.3.3Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

library(lmerTest)    # For hypothesis testing in linear mixed-effects models

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

library(gridExtra)   # For arranging multiple grid graphics

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

library(cowplot)     # For creating complex ggplot2 layouts
library(DescTools)   # For descriptive statistics and data exploration
library(sjmisc)      # For data preparation and variable recoding

Learn more about sjmisc with 'browseVignettes("sjmisc")'.

Attaching package: ‘sjmisc’

The following object is masked from ‘package:DescTools’:

    %nin%

The following object is masked from ‘package:tidyr’:

    replace_na

library(readxl)      # For reading Excel files
library(sjPlot)      # For generating plots and summary tables

Attaching package: ‘sjPlot’

The following objects are masked from ‘package:cowplot’:

    plot_grid, save_plot

library(broom.mixed) # For tidying up model outputs from mixed models
library(scales)      # For scaling and formatting of axes (e.g., percent_format)

Attaching package: ‘scales’

The following object is masked from ‘package:readr’:

    col_factor

library(ggeffects)   # For visualizing effects from regression models

Attaching package: ‘ggeffects’

The following object is masked from ‘package:cowplot’:

    get_title

library(ggridges)    # For creating ridge plots
library(patchwork)   # For combining ggplot2 plots

Attaching package: ‘patchwork’

The following object is masked from ‘package:cowplot’:

    align_plots

library(smplot2)     # For creating summary plots

Warning: package ‘smplot2’ was built under R version 4.3.3Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Updated tutorial for smplot2: smin95.github.io/dataviz/

library(ggpubr)      # For ggplot2-based publication-ready plots

Attaching package: ‘ggpubr’

The following object is masked from ‘package:cowplot’:

    get_legend

library(purrr)       # For functional programming and data manipulation

Attaching package: ‘purrr’

The following object is masked from ‘package:scales’:

    discard

The following object is masked from ‘package:sjmisc’:

    is_empty

The following object is masked from ‘package:metR’:

    cross

1.2. Data Loading

# Read the main data file
df <- read_csv("~/Documents/Research/Code/F1_F0_cue_weighting/Data/data_annotated_v_optim_param_zsco_wth_outliers.csv",show_col_types = FALSE)

# Convert AGE column to factor and specify custom levels
df$AGE <- factor(df$AGE, levels = c("00;06", "00;07", "00;08", "00;09", "00;10", "00;11", "01;00", "01;01", "01;02", "01;03", "01;04", "01;05", "01;06", "01;07", "01;08", "01;09", "01;10", "01;11", "02;00"))


# Read the supplementary data file
cumvoc.y_ALLE_NH <- read_excel("/Users/jeremygenette_studio/Documents/Research/Code/F1_F0_cue_weighting/Data/cumvoc_ALLE_NH_vf.xlsx")

# Process supplementary data
cumvoc.y_ALLE_NH <- cumvoc.y_ALLE_NH %>% 
  filter(!is.na(`Cum voc`)) %>% 
  select(-`Chronage/HearAge`) %>% 
  mutate(item = sub("_.*", "",`OPMERKING: cum op woordvormen, niet op lemma`))

# Merge supplementary data with main data
df <- left_join(df, cumvoc.y_ALLE_NH, by = "item")

# Filter data based on utt_type
df <- df %>%
  filter(utt_type == "LEX") 
# Print the first few rows of the merged and filtered data frame

# Convert pho_vwl_nucl column to factor
df$pho_vwl_nucl <- as.factor(df$pho_vwl_nucl)

# Recode the levels for annotation consistency
df$pho_vwl_nucl <- recode_factor(df$pho_vwl_nucl,
                                  "M" = NA_character_, # diphtong
                                  "L" = NA_character_, # diphtong
                                  "K" = NA_character_, # diphtong
                                  "H" = NA_character_, # glottal stop
                                  "2" = "@", # annotation consistency
                                  "<" = "u", # annotation consistency
                                  ")" = "}") # annotation consistency

# Remove NA levels
df$pho_vwl_nucl <- droplevels(df$pho_vwl_nucl)

# Print the first few rows of the merged and filtered data frame
print(head(df))

2. Data preprocessing ⚙️

2.1. Normalization of F0 based on regression approach from Barreda and Nearey (2018) applied to Lobanov’s (1971) Normalization technique

2.1.1. Data preprocessing

“To implement this analysis in R, it is first assumed that the data are available in a data frame object in a “long” format with only one log-formant measurement per row. Further, it is assumed that each row of the data frame has (at least) four columns, labeled: G for the single formant measurement, V indicating vowel, K indicating formant number, and S indicating the speaker […]an additional variable (N) may be created to represent the Nvk terms, using the R command: N ¼ factor[interaction (V,K)].” (p.507)

# Convert from wide to long_F0 format
df_long_F0 <- df %>%  # Create a new data frame by transforming the existing one
  select(item, CHILD, AGE, ID_vwl, pho_vwl_nucl, F0) %>%  # Select specific columns from the original data frame
  gather(key = "Frequency", value = "Value", F0)  # Reshape the data from wide to long_F0 format, combining F0 and F1 columns into key-value pairs 

df_long_F0$G <- df_long_F0$Value  # assign Value to 'G'
df_long_F0$S <- as.factor(df_long_F0$item)  # Convert the 'item' (=recording session) column to a factor and assign it to 'S'
df_long_F0$V <- as.factor(df_long_F0$pho_vwl_nucl)  # Convert the 'pho_vwl_nucl' column to a factor and assign it to 'V'

2.1.2. Regression

Considering the anticipated variability in how different vowel types affect F0 and F1 frequencies across speakers, which was not addressed in Barreda’s study, we find it imperative to incorporate interactions.

M_F0 = lm(data = df_long_F0,
       formula = G ~ 0 + S * V, contrasts = list(V=contr.sum))

saveRDS(M_F0, "~/Documents/Research/Code/F1_F0_cue_weighting/VF/MF0.rds")

2.1.3. Coefficients extraction

# Extract coefficients for F0 for each speaker
Coefficients_F0 <- summary(M_F0)$coefficients %>%
  as.data.frame() %>%                                          # Convert coefficients to data frame
  filter(grepl("^S", rownames(.))) %>%                         # Filter rows with speaker IDs
  rename(Speaker_Estimate_F0 = Estimate, Speaker_Std_Error_F0 = `Std. Error`) %>%  # Rename columns
  mutate(item = sub("^S", "", rownames(.))) %>%                # Extract speaker IDs
  filter(!grepl(":", item))                                    # Exclude rows with ":" in speaker IDs

# Left join Coefficients_F0 with df_long_F0 by the 'item' column
df_long_with_coef_F0 <- left_join(df_long_F0, Coefficients_F0, by = "item")

# Calculate F0 regression Normalized values
df_long_with_coef_F0 <- df_long_with_coef_F0 %>%
  group_by(item) %>%                                           # Group by speaker IDs
  mutate(F0_Regression_Normalized = (Value - Speaker_Estimate_F0) / (Speaker_Std_Error_F0 * sqrt(n())))  # Calculate Normalized values

# Print head of df_long_with_coef_F0
head(df_long_with_coef_F0)

2.2. Normalization of F1 based on regression approach from Barreda and Nearey (2018) applied to Lobanov’s (1971) Normalization technique

2.2.1. Data preprocessing

“To implement this analysis in R, it is first assumed that the data are available in a data frame object in a “long” format with only one log-formant measurement per row. Further, it is assumed that each row of the data frame has (at least) four columns, labeled: G for the single formant measurement, V indicating vowel, K indicating formant number, and S indicating the speaker […]an additional variable (N) may be created to represent the Nvk terms, using the R command: N ¼ factor[interaction (V,K)].” (p.507)

# Convert from wide to long_F1 format
df_long_F1 <- df %>%  # Create a new data frame by transforming the existing one
  select(item, CHILD, AGE, ID_vwl, pho_vwl_nucl, F1) %>%  # Select specific columns from the original data frame
  gather(key = "Frequency", value = "Value", F1)  # Reshape the data from wide to long_F1 format, combining F1 and F1 columns into key-value pairs 

df_long_F1$G <- df_long_F1$Value  # assign Value to 'G'
df_long_F1$S <- as.factor(df_long_F1$item)  # Convert the 'item' (=recording session) column to a factor and assign it to 'S'
df_long_F1$V <- as.factor(df_long_F1$pho_vwl_nucl)  # Convert the 'pho_vwl_nucl' column to a factor and assign it to 'V'

2.2.2. Regression

Considering the anticipated variability in how different vowel types affect F0 and F1 frequencies across speakers, which was not addressed in Barreda’s study, we find it imperative to incorporate interactions.

M_F1 = lm(data = df_long_F1,
       formula = G ~ 0 + S * V, contrasts = list(V=contr.sum))

saveRDS(M_F1, "~/Documents/Research/Code/F1_F0_cue_weighting/VF/MF1.rds")

2.2.3. Coefficients extraction

# Extract coefficients for F1 for each speaker
Coefficients_F1 <- summary(M_F1)$coefficients %>%
  as.data.frame() %>%                                          # Convert coefficients to data frame
  filter(grepl("^S", rownames(.))) %>%                         # Filter rows with speaker IDs
  rename(Speaker_Estimate_F1 = Estimate, Speaker_Std_Error_F1 = `Std. Error`) %>%  # Rename columns
  mutate(item = sub("^S", "", rownames(.))) %>%                # Extract speaker IDs
  filter(!grepl(":", item))                                    # Exclude rows with ":" in speaker IDs

# Left join Coefficients_F1 with df_long_F1 by the 'item' column
df_long_with_coef_F1 <- left_join(df_long_F1, Coefficients_F1, by = "item")

# Calculate F1 regression Normalized values
df_long_with_coef_F1 <- df_long_with_coef_F1 %>%
  group_by(item) %>%                                           # Group by speaker IDs
  mutate(F1_Regression_Normalized = (Value - Speaker_Estimate_F1) / (Speaker_Std_Error_F1 * sqrt(n())))  # Calculate Normalized values

# Print head of df_long_with_coef_F1
head(df_long_with_coef_F1)

2.3 Gather Normalization of F0 and F1

# Left join df_long_with_coef_F0 with df_long_with_coef_F1 by the 'ID_vwl' column
df_with_coef_F0_F1_lex <- left_join(df_long_with_coef_F0, df_long_with_coef_F1, by = "ID_vwl") %>% 
  rename(item = item.x)  # Rename the 'item' column

# Define the indices of the columns to keep
cols_to_keep <- c(1, 2, 3, 4, 5, 15, 29, 7, 21)

# Extract the desired columns from the dataframe
df_with_coef_F0_F1_lex <- df_with_coef_F0_F1_lex[, cols_to_keep]

# Remove the suffix ".x" from the column names when needed
colnames(df_with_coef_F0_F1_lex) <- gsub("\\.x$", "", colnames(df_with_coef_F0_F1_lex))

df_with_coef_F0_F1_lex$F0 <- df_with_coef_F0_F1_lex$Value
df_with_coef_F0_F1_lex$F1 <- df_with_coef_F0_F1_lex$Value.y

df_with_coef_F0_F1_lex <- df_with_coef_F0_F1_lex %>% 
  select(-c(Value, Value.y))

head(df_with_coef_F0_F1_lex)

2.4. Selection of relevant items

df_combined<- df_with_coef_F0_F1_lex %>% 
  mutate(hgt = case_when(
    pho_vwl_nucl %in% c('A', 'a') ~ 'low',
    pho_vwl_nucl %in% c('u', '<', 'i') ~ 'high',
    TRUE ~ NA_character_  # Handle other cases if needed
  ))  %>% 
  filter(is.na(hgt)==F) %>% 
  select(ID_vwl, F0_Regression_Normalized, F1_Regression_Normalized, pho_vwl_nucl, hgt, CHILD, item, F0, F1)

2.5. Preparation df for analysis

# Converting 'hgt' to a factor
df_combined$hgt <- as.factor(df_combined$hgt)

# Categorizing 'frt' based on 'pho_vwl_nucl'
df_combined$frt[df_combined$pho_vwl_nucl %in% c("A", "u", "<")] <- "bck"

Warning: Unknown or uninitialised column: `frt`.

df_combined$frt[df_combined$pho_vwl_nucl %in% c("i", "a")] <- "frt"

# Merging supplementary data with main data and filtering based on a condition
df_combined <- left_join(df_combined, cumvoc.y_ALLE_NH, by = "item") 
df_combined <- df_combined %>% 
  filter(`Cum voc` > 0)

# Calculating the logarithm of 'cumvoc.y'
df_combined$log_cum <- log(df_combined$`Cum voc`)

# Final dataframe for analysis
df <- df_combined %>% 
  mutate(hgt = as.factor(hgt),
         CHILD = as.factor(CHILD),
         item = as.factor(item),
         frt = as.factor(frt))

 
# Apply deviation (sum) contrasts to 'hgt' and 'frt'
# This codes each level relative to the overall mean of the dependent variable
df$hgt <- factor(df$hgt, levels = c("low", "high"))  # Ensure 'hgt' is a factor with specified levels
contrasts(df$hgt) <- contr.sum(levels(df$hgt))  # Apply contr.sum to 'hgt'
contrasts(df$frt) <- contr.sum(levels(df$frt))  # Apply contr.sum to 'frt'

3. Data analysis 🧮

3.1. Descriptive Statistics

# Calculate mean and standard deviation of F1 for each pho_vwl_nucl category (Normalized)
df_descr <- df %>% 
  group_by(hgt) %>%   # Group by pho_vwl_nucl category
  mutate(F1_F0_Hz = F1-F0,
         F1_F0_Normalized = F1_Regression_Normalized-F0_Regression_Normalized) %>% 
  summarise(
    m_F1_Hz = round(mean(F1), 2),             # Mean F1 (Hz)
    sd_F1_Hz = round(sd(F1), 2),               # Standard deviation of F1 (Hz)
    m_F0_Hz = round(mean(F0), 2),                # Mean F0 (Hz)
    sd_F0_Hz = round(sd(F0), 2),                 # Standard deviation of F0 (Hz)
    m_F1_F0_Hz = round(mean(F1_F0_Hz), 2),          # Mean F1-F0 difference (Hz)
    sd_F1_F0_Hz = round(sd(F1_F0_Hz), 2),           # Standard deviation of F1-F0 difference (Hz)
    m_F1_Normalized = round(mean(F1_Regression_Normalized), 2),     # Mean F1 (Normalized)
    sd_F1_Normalized = round(sd(F1_Regression_Normalized), 2),       # Standard deviation of F1 (Normalized)
    m_F0_Normalized = round(mean(F0_Regression_Normalized), 2),       # Mean F0 (Normalized)
    sd_F0_Normalized = round(sd(F0_Regression_Normalized), 2),       # Standard deviation of F0 (Normalized)
    m_F1_F0_Normalized = round(mean(F1_F0_Normalized), 2),           # Mean F1-F0 difference (Normalized)
    sd_F1_F0_Normalized = round(sd(F1_F0_Normalized), 2)) %>%       # Standard deviation of F1-F0 difference (Normalized)
  filter(!is.na(hgt))   # Remove rows where hgt is NA

print(df_descr)

df_descr_phon <- df %>% 
  group_by( CHILD, pho_vwl_nucl) %>% 
  summarise(counts = n())

`summarise()` has grouped output by 'CHILD'. You can override using the `.groups` argument.

print(df_descr_phon)

df_descr_phon_avg <- df %>% 
  group_by( CHILD, pho_vwl_nucl) %>% 
  summarise(counts = n()) %>% 
  group_by(pho_vwl_nucl) %>%
  summarise(m_counts=mean(counts),
           sd=sd(counts))

`summarise()` has grouped output by 'CHILD'. You can override using the `.groups` argument.

print(df_descr_phon_avg)

NA

3.1.2. Graphs experimental data

# vowel distribution per child plot
df_descr_phon$pho_vwl_nucl <- factor(
  df_descr_phon$pho_vwl_nucl,
  levels = c("i", "u", "a", "A")
)

df_descr_phon <- df_descr_phon %>%
  ungroup() %>%   
  mutate(
    CHILD_anon = factor(
      CHILD,
      labels = paste0("ID",  seq_along(levels(CHILD)))
    )
  )

child_key <- df_descr_phon %>%
  ungroup() %>%
  distinct(CHILD) %>%
  arrange(CHILD) %>%   # optional: alphabetic order
  mutate(
    CHILD_anon = paste0("ID",  seq_along(levels(CHILD)))
  )

count <- ggplot(df_descr_phon,
       aes(x = CHILD_anon, y = counts,
           fill = pho_vwl_nucl, colour = pho_vwl_nucl)) +
  geom_col(alpha = 0.85) +
  labs(
    x = "Child",
    y = "Count",
    fill = "Vowel"
  ) +
  scale_fill_manual(
    values = c(
  "i" = "#1F77B4",  # blue
  "u" = "#17BECF",  # cyan
  "a" = "#FF7F0E",  # orange
  "A" = "#E41A1C"   # red
    ),
    labels = c(
      "i" = "[i]",
      "u" = "[u]",
      "a" = "[a]",
      "A" = "[ɑ]"
    )
  ) +
  scale_colour_manual(
    values = c(
  "i" = "#1F77B4",  # blue
  "u" = "#17BECF",  # cyan
  "a" = "#FF7F0E",  # orange
  "A" = "#E41A1C"   # red
    ),
    guide = "none"   
  ) +
  theme_minimal() +
  theme(
    plot.title = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_line(color = "grey85"),
axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_text(size = 14),
    legend.text  = element_text(size = 16)
  )

count 
ggsave(filename = "~/Desktop/image_eps.eps", plot = count,
       width = 15, height = 6, dpi = 1000)

# F1 density plot
p1 <- ggplot(data = df) +
  geom_density(aes(x = F1_Regression_Normalized, fill = hgt), alpha = 0.5) +
  scale_fill_manual(values = c("low" = "#9F3400", "high" = "#FFBE9F"), name = "Height") +
  labs(x = "F1 Normalized", y = "Density") +
  scale_x_continuous(limits = c(-1.2, 1.2)) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank()  ,  legend.text = element_text(size = 14)      # Increase legend title size (optional)
  )
  

# F0 density plot
p2 <- ggplot(data = df) +
  geom_density(aes(x = F0_Regression_Normalized, fill = hgt), alpha = 0.5) +
  scale_fill_manual(values = c("low" = "#570987", "high" = "#D397F8"), name = "Height") +
  labs(x = "F0 Normalized", y = "Density") +
  scale_x_continuous(limits = c(-1.2, 1.2)) +
  theme_minimal() +
  theme(
    legend.position = "bottom",  # Position the legend on the right side,
    panel.grid.major.y = element_blank(),
   
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),  legend.text = element_text(size = 14)  
  )

# Combine with shared legend
combined_plots <- (p1 + p2) +
 plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

combined_plots

ggsave(filename = "~/Desktop/image_eps.eps", plot = combined_plots,
       width = 10, height = 6, dpi = 1000)

df$pho_vwl_nucl <- factor(
  df$pho_vwl_nucl,
  levels = c("i", "u", "a", "A")
)

f0_raw <- ggplot(df, aes(x = log_cum, y = F0_Regression_Normalized, group = pho_vwl_nucl)) +
  
  # Ribbon layer (legend shows only fill)
  geom_smooth(
    aes(fill = pho_vwl_nucl),
    method = "loess",
    linewidth = 1.2,
    color = NA,            # no line in this layer
    show.legend = TRUE
  ) +
  
  # Line layer (no legend)
  geom_smooth(
    aes(color = pho_vwl_nucl),
    method = "loess",
    linewidth = 1.2,
    fill = NA,             # don't fill this layer
    show.legend = FALSE
  ) +
  
  # Fill scale controls the legend
  scale_fill_manual(
    name = "Vowel",        # <- legend title
    values = c(
  "i" = "#1F77B4",  # blue
  "u" = "#17BECF",  # cyan
  "a" = "#FF7F0E",  # orange
  "A" = "#E41A1C"   # red
    ),
    labels = c(
      "i" = "[i]",
      "u" = "[u]",
      "a" = "[a]",
      "A" = "[ɑ]"
    )
  ) +
  
  scale_color_manual(
    values = c(
  "i" = "#1F77B4",  # blue
  "u" = "#17BECF",  # cyan
  "a" = "#FF7F0E",  # orange
  "A" = "#E41A1C"   # red
    ),
    guide = "none"           # hide the color legend
  ) +
  
  # Axis labels
  labs(
    x = "log(cumulative vocabulary)",
    y = "Normalized F0"
  ) +
  
  ylim(-1, 1) +
  theme_cowplot() +
  theme(
    legend.position = "right",
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 16),  # bigger legend title
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  )

f0_raw
ggsave(filename = "~/Desktop/f0_raw_eps.eps", plot = f0_raw,
       width = 6, height = 6, dpi = 1000)

df$pho_vwl_nucl <- factor(
  df$pho_vwl_nucl,
  levels = c("i", "u", "a", "A")
)

f1_raw <- ggplot(df, aes(x = log_cum, y = F1_Regression_Normalized, group = pho_vwl_nucl)) +
  
  # Ribbon layer (legend shows only fill)
  geom_smooth(
    aes(fill = pho_vwl_nucl),
    method = "loess",
    linewidth = 1.2,
    color = NA,            # no line in this layer
    show.legend = TRUE
  ) +
  
  # Line layer (no legend)
  geom_smooth(
    aes(color = pho_vwl_nucl),
    method = "loess",
    linewidth = 1.2,
    fill = NA,             # don't fill this layer
    show.legend = FALSE
  ) +
  
  # Fill scale controls the legend
  scale_fill_manual(
    name = "Vowel",        # <- legend title
    values = c(
  "i" = "#1F77B4",  # blue
  "u" = "#17BECF",  # cyan
  "a" = "#FF7F0E",  # orange
  "A" = "#E41A1C"   # red
    ),
    labels = c(
      "i" = "[i]",
      "u" = "[u]",
      "a" = "[a]",
      "A" = "[ɑ]"
    )
  ) +
  
  scale_color_manual(
    values = c(
  "i" = "#1F77B4",  # blue
  "u" = "#17BECF",  # cyan
  "a" = "#FF7F0E",  # orange
  "A" = "#E41A1C"   # red
    ),
    guide = "none"           # hide the color legend
  ) +
  
  # Axis labels
  labs(
    x = "log(cumulative vocabulary)",
    y = "Normalized F1"
  ) +
  
  ylim(-1, 1) +
  theme_cowplot() +
  theme(
    legend.position = "right",
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 16),  # bigger legend title
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  )

f1_raw

ggsave(filename = "~/Desktop/f1_raw.eps", plot = f1_raw,
       width = 6, height = 6, dpi = 1000)

3.2. MLM Modelling

3.2.1. Model 1: Base model

m1 <- glmer(data= df, 
            formula = hgt ~ 1+ 
              (1|CHILD),  family= "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m1, file = "m1.rds")

3.2.2. Model 2: Fixed effect of cumulative vocabulary

m2 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              (1|CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m2, file = "m2.rds")

as.data.frame(anova(m1,m2))

➡️ Better datafit with m2

3.2.3. Model 3: Random effect of cumulative vocabulary

m3 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              (1 + log_cum | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m3, file = "m3.rds")

as.data.frame(anova(m2,m3))

➡️ Better datafit with m3

3.2.4. Model 4: Fixed effect of F1

m4 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              (1 + log_cum | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m4, file = "m4.rds")

as.data.frame(anova(m3,m4))

➡️ Better datafit with m4

3.2.5. Model 5: Random slope effect of F1

m5 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              (1 + log_cum + F1_Regression_Normalized | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m5, file = "m5.rds")

as.data.frame(anova(m4,m5))

➡️ Better datafit with m5

3.2.6. Model 6: Fixed effect of F0

m6 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              (1 + log_cum + F1_Regression_Normalized | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m6, file = "m6.rds")

as.data.frame(anova(m5,m6))

➡️ Better datafit with m6

3.2.7. Model 7: Random slope effect of F0

m7 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m7, file = "m7.rds")

as.data.frame(anova(m6,m7))

➡️ Better datafit with m7

3.2.8. Model 8: Fixed effect of place of articulation

m8 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m8, file = "m8.rds")

as.data.frame(anova(m7,m8))

➡️ Better datafit with m8

3.2.9. Model 9: Random slope effect of place of articulation

m9 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m9, file = "m9.rds")

as.data.frame(anova(m8,m9))

➡️ Better datafit with m9

3.2.10. Model 10: Fixed quadratic effect of cumulative vocabulary interaction

m10 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              I(log_cum^2)+ 
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m10, file = "m10.rds")

as.data.frame(anova(m9,m10))

❌ No better datafit with m10

3.2.11. Model 11: Fixed effect of F0-F1 interaction

m11 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              F1_Regression_Normalized : F0_Regression_Normalized+
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m11, file = "m11.rds")


as.data.frame(anova(m9,m11))

➡️ Better datafit with m11

3.2.12. Model 12: Fixed effect of place of articulation-F1 interaction

m12 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              F1_Regression_Normalized : F0_Regression_Normalized+
              F1_Regression_Normalized : frt +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m12, file = "m12.rds")

as.data.frame(anova(m11,m12))

➡️ Better datafit with m12

3.2.13. Model 13: Fixed effect of cumulative vocabulary-F1 interaction

m13 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              F1_Regression_Normalized : F0_Regression_Normalized+
              F1_Regression_Normalized : frt +
              F1_Regression_Normalized : log_cum +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m13, file = "m13.rds")

as.data.frame(anova(m12,m13))

➡️ Better with m13

3.2.14. Model 14: Fixed effect of F0-place of articulation interaction

m14 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              F1_Regression_Normalized : F0_Regression_Normalized+
              F1_Regression_Normalized : frt +
              F1_Regression_Normalized : log_cum +
              F0_Regression_Normalized : frt +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m14, file = "m14.rds")

as.data.frame(anova(m13,m14))

➡️ Better datafit with m14

3.2.15. Model 15: Fixed effect of cumulative vocabulary-F0 interaction

m15 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              F1_Regression_Normalized : F0_Regression_Normalized+
              F1_Regression_Normalized : frt +
              F1_Regression_Normalized : log_cum +
              F0_Regression_Normalized : frt +
              F0_Regression_Normalized : log_cum +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m15, file = "m15.rds")

as.data.frame(anova(m14,m15))

➡️ Better datafit with m15

3.2.16. Model 16: Fixed effect of cumulative vocabulary-F0-F1 interaction

m16 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              F1_Regression_Normalized : F0_Regression_Normalized+
              F1_Regression_Normalized : frt +
              F1_Regression_Normalized : log_cum +
              F0_Regression_Normalized : frt +
              F0_Regression_Normalized : log_cum +
              F1_Regression_Normalized : F0_Regression_Normalized : log_cum +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m16, file = "m16.rds")

as.data.frame(anova(m15,m16))

❌ No better datafit with m16

3.2.17. Model 17: Fixed effect of place of articulation-F0-F1 interaction

m17 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              F1_Regression_Normalized : F0_Regression_Normalized+
              F1_Regression_Normalized : frt +
              F1_Regression_Normalized : log_cum +
              F0_Regression_Normalized : frt +
              F0_Regression_Normalized : log_cum +
              F1_Regression_Normalized : F0_Regression_Normalized : frt +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m17, file = "m17.rds")

as.data.frame(anova(m15,m17))

➡️ Better datafit with m17

3.2.18. Model 18: Fixed effect of cumulative vocabulary-place of articulation-F0-F1 interaction

m18 <- glmer(data = df, 
            formula = hgt ~ log_cum + 
              F1_Regression_Normalized+
              F0_Regression_Normalized+
              frt +
              F1_Regression_Normalized : F0_Regression_Normalized+
              F1_Regression_Normalized : frt +
              F1_Regression_Normalized : log_cum +
              F0_Regression_Normalized : frt +
              F0_Regression_Normalized : log_cum +
              F1_Regression_Normalized : F0_Regression_Normalized : frt +
              F1_Regression_Normalized : F0_Regression_Normalized : frt : log_cum +
              (1 + log_cum + F1_Regression_Normalized + F0_Regression_Normalized + frt | CHILD),  family = "binomial"(link = "logit"),
            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2000000)))
saveRDS(m18, file = "m18.rds")

as.data.frame(anova(m17,m18))

❌ No better datafit with m18

3.3. Final Model ✅

final_model <- m17
print(summary(final_model))

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: hgt ~ log_cum + F1_Regression_Normalized + F0_Regression_Normalized +  
    frt + F1_Regression_Normalized:F0_Regression_Normalized +  
    F1_Regression_Normalized:frt + F1_Regression_Normalized:log_cum +  
    F0_Regression_Normalized:frt + F0_Regression_Normalized:log_cum +  
    F1_Regression_Normalized:F0_Regression_Normalized:frt + (1 +  
    log_cum + F1_Regression_Normalized + F0_Regression_Normalized +      frt | CHILD)
   Data: df
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+06))

     AIC      BIC   logLik deviance df.resid 
  9029.2   9223.3  -4488.6   8977.2    12880 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
 -6.281  -0.375  -0.156   0.248 240.387 

Random effects:
 Groups Name                     Variance Std.Dev. Corr                   
 CHILD  (Intercept)              2.2732   1.5077                          
        log_cum                  0.1064   0.3263   -0.97                  
        F1_Regression_Normalized 2.8553   1.6898    0.33 -0.44            
        F0_Regression_Normalized 3.0895   1.7577   -0.01  0.05 -0.02      
        frt1                     0.1505   0.3880    0.28 -0.32  0.34 -0.15
Number of obs: 12906, groups:  CHILD, 30

Fixed effects:
                                                       Estimate Std. Error z value Pr(>|z|)
(Intercept)                                            -2.00898    0.32324  -6.215 5.13e-10
log_cum                                                 0.27316    0.07079   3.859 0.000114
F1_Regression_Normalized                               -1.87675    0.70340  -2.668 0.007628
F0_Regression_Normalized                                1.66602    0.59373   2.806 0.005015
frt1                                                    0.52375    0.07874   6.651 2.91e-11
F1_Regression_Normalized:F0_Regression_Normalized      -1.87024    0.62241  -3.005 0.002657
F1_Regression_Normalized:frt1                           0.79509    0.17455   4.555 5.23e-06
log_cum:F1_Regression_Normalized                       -1.42642    0.14866  -9.595  < 2e-16
F0_Regression_Normalized:frt1                          -0.94897    0.14672  -6.468 9.95e-11
log_cum:F0_Regression_Normalized                        0.40528    0.11618   3.488 0.000486
F1_Regression_Normalized:F0_Regression_Normalized:frt1  1.16619    0.58946   1.978 0.047883
                                                          
(Intercept)                                            ***
log_cum                                                ***
F1_Regression_Normalized                               ** 
F0_Regression_Normalized                               ** 
frt1                                                   ***
F1_Regression_Normalized:F0_Regression_Normalized      ** 
F1_Regression_Normalized:frt1                          ***
log_cum:F1_Regression_Normalized                       ***
F0_Regression_Normalized:frt1                          ***
log_cum:F0_Regression_Normalized                       ***
F1_Regression_Normalized:F0_Regression_Normalized:frt1 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
               (Intr) log_cm F1_Rg_N F0_Rg_N frt1   F1_Rg_N:F0_R_N F1_R_N:1 l_:F1_ F0_R_N:
log_cum        -0.970                                                                     
F1_Rgrssn_N     0.084 -0.122                                                              
F0_Rgrssn_N    -0.002  0.017 -0.125                                                       
frt1            0.202 -0.234  0.158  -0.079                                               
F1_Rg_N:F0_R_N -0.056  0.047  0.109  -0.181   0.024                                       
F1_Rgrs_N:1     0.006 -0.020 -0.041  -0.002  -0.010  0.035                                
lg_c:F1_R_N     0.054 -0.057 -0.855   0.134  -0.039 -0.184          0.014                 
F0_Rgrs_N:1     0.021 -0.017  0.010  -0.049  -0.111  0.000         -0.141    0.004        
lg_c:F0_R_N    -0.011 -0.001  0.128  -0.786   0.015  0.176          0.018   -0.158  0.032 
F1_R_N:F0_R_N:  0.007  0.000 -0.012   0.135  -0.037 -0.197         -0.245    0.012 -0.203 
               l_:F0_
log_cum              
F1_Rgrssn_N          
F0_Rgrssn_N          
frt1                 
F1_Rg_N:F0_R_N       
F1_Rgrs_N:1          
lg_c:F1_R_N          
F0_Rgrs_N:1          
lg_c:F0_R_N          
F1_R_N:F0_R_N: -0.171

# Assuming final_model is your mixed-effects model
summary_final_model <- summary(final_model)

formula <- final_model@call$formula

# Fixed effects summary
print(as.data.frame(round(coef(summary_final_model),3)))

print(as.data.frame(coef(summary_final_model)))

4. Data vizaulization 📊

4.0. Refit model with cum.voc. centred at other values

# Center by subtracting the mean of 'log_cum'
df_no_ctr <- df
df_no_ctr$log_cum <- df_no_ctr$log_cum - mean(df_no_ctr$log_cum)  # Center by the mean

# Fit the model again using the centered log_cum (mean-centered)
refit_model_ctr_mean <- glmer(
  formula = final_model@call$formula,
  data = df_no_ctr,
  family = binomial,
  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
)

df_max <- df
df_max$log_cum <- df_max$log_cum-max(df_max$log_cum)

refit_model_max <- glmer(
formula <- final_model@call$formula,
  data = df_max,
  family = binomial,
  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
)

4.1. F0 and F1 effects on vowel height

# Plot F1 with separate lines for each level of facto_hgt
plot_f1 <- plot_model(final_model, 
                      type = "eff", 
                      terms = c("F1_Regression_Normalized [all]", "frt"),  
                      mdrt.values = "meansd", 
                      color = "orange") +  # color lines by factor
  theme_cowplot() +
  xlim(c(-2, 2)) +
  labs(y = "Probability of a vowel being high", 
       x = "Normalized F1") +
  theme(plot.title = element_blank(), 
        plot.subtitle = element_blank(), 
        plot.caption = element_blank())+
  aes(linetype = group_col)+
  scale_linetype_discrete(name = "Place of \narticulation",labels = c("back", "front"))  +
  scale_color_manual(values = c("orange","orange"), guide = "none")+
  scale_fill_manual(values = c("orange","orange"), guide = "none")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

plot_f0 <- plot_model(final_model, 
                      type = "eff", 
                      terms = c("F0_Regression_Normalized [all]", "frt"), 
                      mdrt.values = "meansd", 
                      at = list(log_cum = max(log_cum)),
                      color = "purple") +
    aes(linetype = group_col)+
  theme_cowplot() +
  xlim(c(-2, 2)) +
  labs(y = "", x = "Normalized F0") +
  theme(plot.title = element_blank(), 
        plot.subtitle = element_blank(), 
        plot.caption = element_blank())+
  scale_linetype_discrete(name = "Place of \narticulation",  labels = c("back", "front"))  +
  scale_color_manual(values = c("purple","purple"), guide = "none")+
  scale_fill_manual(values = c("purple","purple"), guide = "none")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

# Combine plots
combined_plots <- grid.arrange(plot_f1, plot_f0, ncol = 2)

# Print combined plot
print(combined_plots)

TableGrob (1 x 2) "arrange": 2 grobs

# Save
ggsave(filename = "~/Desktop/f1_vs_f0_by_frt_eps.eps", 
       plot = combined_plots, 
       width = 10, height = 6, dpi = 1000)

Warning: semi-transparency is not supported on this device: reported only once per page

4.2. Developmental trajectories

# Define colors for F1 and F0
f1_color <- c("#9F3400","#FF681F","#FFBE9F")
f0_color <- c("#570987","#A020F0","#D397F8")

# Plot for F1
plot_f1 <- plot_model(final_model, type = "eff", terms = c("log_cum [all]", "F1_Regression_Normalized")) +
  labs(color = "Normalized F1", y = "Predicted probability of a vowel being high", x = "log(cumulative vocabulary)", title = "") +
  theme_cowplot() +
  scale_color_manual(values = f1_color) +  # Specify color for F1
  scale_fill_manual(values = f1_color) +  # Specify color for F1
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +    # Set y-axis limits from 0 to 1
  theme(legend.position = "bottom", legend.box = "horizontal")  

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

# Plot for F0
plot_f0 <- plot_model(final_model, type = "eff", terms = c("log_cum [all]", "F0_Regression_Normalized")) +
  labs(color = "Normalized F0", y = "Predicted probability of a vowel being high", x = "log(cumulative vocabulary)", title = "") +
  theme_cowplot() +
  scale_color_manual(values = f0_color) +  # Specify color for F0
  scale_fill_manual(values = f0_color) +  # Specify color for F0
   scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +    # Set y-axis limits from 0 to 1
  theme(legend.position = "bottom", legend.box = "horizontal")  

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

# Adjust plot width and height
plot_f1 

plot_f0


ggsave(
  filename = "~/Desktop/plot_f1_eps.eps",
  plot = plot_f1,
  width = 6,
  height = 5,
  units = "in"
)

ggsave(
  filename = "~/Desktop/plot_f0_eps.eps",
  plot = plot_f0,
  width = 6,
  height = 5,
  units = "in"
)

df_plot_f1 <- plot_f1$data
df_plot_f0 <- plot_f0$data

print(as.data.frame(df_plot_f1))

print(as.data.frame(df_plot_f0))

pr_F1<-predict_response(final_model,c("F1_Regression_Normalized", "log_cum" ))

Data were 'prettified'. Consider using `terms="F1_Regression_Normalized [all]"` to
  get smooth plots.

jn_pr_F1 <- plot(johnson_neyman(pr_F1))  +
    scale_color_manual(values = c("inconsistent" = "grey", "positive/negative" = "orange")) +
  scale_fill_manual(values = c("inconsistent" = "grey", "positive/negative" = "orange")) +
  labs(y = "Slope of Normalized F1", x = "log(cumulative vocabulary)", title = "") +
  theme_cowplot() +
  theme(legend.position = "bottom", legend.box = "horizontal") + 
  ylim(c(-0.5,0.5))

Warning: For this model type, `marginaleffects` only takes into account the uncertainty in
  fixed-effect parameters. You can use the `re.form=NA` argument to acknowledge this
  explicitly and silence this warning.Warning: For this model type, `marginaleffects` only takes into account the uncertainty in
  fixed-effect parameters. You can use the `re.form=NA` argument to acknowledge this
  explicitly and silence this warning.

The association between `F1_Regression_Normalized` and `hgt` is negative for values
  of `log_cum` higher than 0.32. There were no clear associations for values of
  `log_cum` lower than 0.32.

pr_F0<-predict_response(final_model,c("F0_Regression_Normalized", "log_cum"))

Data were 'prettified'. Consider using `terms="F0_Regression_Normalized [all]"` to
  get smooth plots.

jn_pr_F0 <- plot(johnson_neyman(pr_F0))  +
  scale_color_manual(values = c("inconsistent" = "grey", "positive/negative" = "purple")) +
  scale_fill_manual(values = c("inconsistent" = "grey", "positive/negative" = "purple")) +
  labs(y = "Slope of Normalized F0", x = "log(cumulative vocabulary)", title = "") +
  theme_cowplot() +
  theme(legend.position = "bottom", 
        legend.box = "horizontal")+ 
  ylim(c(-0.5,0.5))

Warning: For this model type, `marginaleffects` only takes into account the uncertainty in
  fixed-effect parameters. You can use the `re.form=NA` argument to acknowledge this
  explicitly and silence this warning.Warning: For this model type, `marginaleffects` only takes into account the uncertainty in
  fixed-effect parameters. You can use the `re.form=NA` argument to acknowledge this
  explicitly and silence this warning.

The association between `F0_Regression_Normalized` and `hgt` is positive for values
  of `log_cum` higher than 0.85. There were no clear associations for values of
  `log_cum` lower than 0.85.

jn_pr_F1

jn_pr_F0

ggsave(filename = "~/Desktop/jn_pr_F0_eps.eps", 
       plot = jn_pr_F0, 
       width = 10, height = 6, dpi = 1000)

ggsave(filename = "~/Desktop/jn_pr_F1_eps.eps", 
       plot = jn_pr_F1, 
       width = 10, height = 6, dpi = 1000)

4.3. Correlation between cues

corr_full <- cor.test(df$F1_Regression_Normalized,
                      df$F0_Regression_Normalized, method = 'pearson')

high <- df %>% filter(hgt == "high" )
low  <- df %>% filter(hgt == "low")

corr_high <- cor.test(high$F1_Regression_Normalized,
                      high$F0_Regression_Normalized)

corr_low <- cor.test(low$F1_Regression_Normalized,
                     low$F0_Regression_Normalized)

# Put all results in a list
corr_results <- list(
  full = corr_full,
  high = corr_high,
  low  = corr_low
)

corr_results

$full

    Pearson's product-moment correlation

data:  df$F1_Regression_Normalized and df$F0_Regression_Normalized
t = 10.484, df = 12904, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.0747654 0.1089797
sample estimates:
       cor 
0.09189966 


$high

    Pearson's product-moment correlation

data:  high$F1_Regression_Normalized and high$F0_Regression_Normalized
t = 5.6985, df = 3349, p-value = 1.313e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.06434937 0.13142014
sample estimates:
       cor 
0.09799603 


$low

    Pearson's product-moment correlation

data:  low$F1_Regression_Normalized and low$F0_Regression_Normalized
t = 21.73, df = 9553, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1978370 0.2360515
sample estimates:
      cor 
0.2170274 

df_avg <- df %>%
  group_by(CHILD,hgt) %>%
  summarise(
    F1_avg = mean(F1_Regression_Normalized, na.rm = TRUE),
    F0_avg = mean(F0_Regression_Normalized, na.rm = TRUE),
    .groups = "drop"
  )

corr_full_avg <- cor.test(df_avg$F1_avg, df_avg$F0_avg, method = "pearson")

df_avg <- df %>%
  group_by(CHILD,hgt) %>%
  summarise(
    F1_avg = mean(F1_Regression_Normalized, na.rm = TRUE),
    F0_avg = mean(F0_Regression_Normalized, na.rm = TRUE),
    .groups = "drop"
  )

high_avg <- df_avg %>% filter(hgt == "high")
low_avg  <- df_avg %>% filter(hgt == "low")

corr_high_avg <- cor.test(high_avg$F1_avg, high_avg$F0_avg, method = "pearson")
corr_low_avg  <- cor.test(low_avg$F1_avg,  low_avg$F0_avg,  method = "pearson")


corr_results_avg <- list(
  full_avg = corr_full_avg,
  high_avg = corr_high_avg,
  low_avg  = corr_low_avg
)

corr_results_avg

$full_avg

    Pearson's product-moment correlation

data:  df_avg$F1_avg and df_avg$F0_avg
t = -4.3161, df = 58, p-value = 6.269e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6638678 -0.2733592
sample estimates:
       cor 
-0.4930594 


$high_avg

    Pearson's product-moment correlation

data:  high_avg$F1_avg and high_avg$F0_avg
t = -0.16062, df = 28, p-value = 0.8735
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.3863857  0.3335755
sample estimates:
        cor 
-0.03033985 


$low_avg

    Pearson's product-moment correlation

data:  low_avg$F1_avg and low_avg$F0_avg
t = 1.4681, df = 28, p-value = 0.1532
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1028207  0.5724807
sample estimates:
      cor 
0.2673521 

# Extract coefficient values and convert to a data frame for mean-centered model
coefficients_mean <- as.data.frame(coef(refit_model_ctr_mean )$CHILD)
coefficients_mean$id <- row.names(coefficients_mean)
coefficients_mean$model <- "mean_centered"

# Combine final, refit, and mean-centered coefficients into one data frame
coefficients_combined <- rbind(coefficients_mean)

# Make sure id is a factor
coefficients_combined$id <- as.factor(coefficients_combined$id)

# Gather F1 and F0 into one 'variable' column
coefficients_long <- coefficients_combined %>%
  select(id, model, F1_Regression_Normalized, F0_Regression_Normalized) %>%
  gather(key = "variable", value = "value", F1_Regression_Normalized, F0_Regression_Normalized)

# Spread by model (final, refit, and mean-centered)
coefficients_wide <- coefficients_long %>%
  spread(key = model, value = value)

# Pivot wider to include final, refit, and mean-centered values
coefficients_final_wide <- coefficients_wide %>%
  pivot_wider(
    names_from = variable,
    values_from = c(mean_centered)
  )

plot <- ggplot(coefficients_final_wide, aes(x =F1_Regression_Normalized,
                                    y = F0_Regression_Normalized)) +
  geom_point(color = "grey60", size = 3, alpha = 0.8) +
  sm_statCorr(color = "#5F021F")+
  xlim(-12,0)+
  ylim(0,12)+
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "horizontal",
    panel.grid.major = element_line(color = "gray95", size = 0.5),  # Subtle grid lines
    panel.grid.minor = element_blank(),  # Remove minor grid lines for a cleaner look
    axis.title = element_text(size = 14, face = "bold"),
    axis.text = element_text(size = 12),
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),,
    plot.margin = margin(20, 20, 20, 20)  # Add margin for aesthetics
  ) +
  labs(
    x = "F1 Regression Normalized",
    y = "F0 Regression Normalized")+
  coord_fixed()

plot

ggsave(filename = "~/Desktop/corr_eps.eps", 
       plot = plot, 
       width = 10, height = 6, dpi = 1000)

4.4. Between children variability

# Extract coefficient values and convert to a data frame
coefficients_0<- as.data.frame(coef(final_model)$CHILD)
coefficients_0$id <- row.names(coefficients_0)   # Add ID as a column
coefficients_0$id <- as.factor(coefficients_0$id)   # Convert ID to factor

# Reshape the data frame from wide to long format
coefficients_0<- pivot_longer(coefficients_0, cols = c(F0_Regression_Normalized, F1_Regression_Normalized), names_to = "Variable", values_to = "Value")

# Extract coefficient values and convert to a data frame
coefficients_max<- as.data.frame(coef(refit_model_max)$CHILD)
coefficients_max$id <- row.names(coefficients_max)   # Add ID as a column
coefficients_max$id <- as.factor(coefficients_max$id)   # Convert ID to factor


# Reshape the data frame from wide to long format
coefficients_max<- pivot_longer(coefficients_max, cols = c(F0_Regression_Normalized, F1_Regression_Normalized), names_to = "Variable", values_to = "Value")

# Add a "Centering" column to identify where the coefficients come from
coefficients_0$Centering <- "Centered_at_0"
coefficients_max$Centering <- "Centered_at_Max"


# Combine the two data frames
coef_full <- bind_rows(coefficients_0, coefficients_max)

# Add a new column to control alpha
coef_full$Alpha <- ifelse(coef_full$Centering == "Centered_at_0", 0.9, 1)  # Circles = 0.3, Triangles = 1

# Update the levels and labels for better readability
coef_full$Variable <- factor(coef_full$Variable, 
                             levels = c("F0_Regression_Normalized", "F1_Regression_Normalized"),
                             labels = c("F0", "F1"))

coef_full$Centering <- factor(coef_full$Centering,
                              levels = c("Centered_at_0", "Centered_at_Max"),
                              labels = c("At cumulative vocabulary = 0", "At maximum cumulative vocabulary"))


# Reorder ids based on F1 values (choose which centering you want)
coef_full_ordered <- coef_full %>%
  filter(Variable == "F0") %>%
  arrange(Value) %>%
  mutate(id = factor(id, levels = unique(id))) %>%
  select(id) %>%
  right_join(coef_full, by = "id") %>%   # restore full dataset with ordering
  left_join(
    df_descr_phon %>% distinct(id = CHILD, label = CHILD_anon),
    by = "id"
  )

Warning: Detected an unexpected many-to-many relationship between `x` and `y`.

coef_full_ordered <- coef_full_ordered %>%
  mutate(label = factor(label, levels = unique(label[order(id)])))

plot <- ggplot(coef_full_ordered, aes(x = label, y = Value)) +
  geom_point(size = 3, aes(color = Variable, shape = Centering)) +
  scale_shape_manual(values = c(
    "At cumulative vocabulary = 0" = 1,
    "At maximum cumulative vocabulary" = 16
  )) +
  scale_color_manual(values = c("F0" = "purple", "F1" = "orange")) +
  geom_segment(aes(x = label, xend = label, y = 0, yend = Value, color = Variable),
               size = 0.5) +
  labs(
    title = "",
    x = "Child label",
    y = "Individual Coefficient Value",
    color = "Cue",
    shape = "Cumulative vocabulary"
  ) +
  #scale_x_discrete(labels = label) +
  theme_minimal() +
  theme(
    legend.position = "right",
    legend.direction = "vertical",
    legend.box = "vertical",
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 16, hjust = 0.5),
    panel.grid.minor.y = element_blank(),
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 13)
  )

plot

ggsave(filename = "~/Desktop/lolliplot_eps.eps", 
       plot = plot, 
       width = 20, height = 6, dpi = 1000)

# Choose constant values for F1_Regression_Normalized and F0_Regression_Normalized
constant_F1_low <- mean(df$F1_Regression_Normalized) - 1*sd(df$F1_Regression_Normalized)
constant_F0_high <- mean(df$F0_Regression_Normalized) + 1*sd(df$F0_Regression_Normalized)
constant_F1_high <- mean(df$F1_Regression_Normalized) + 1*sd(df$F1_Regression_Normalized)
constant_F0_low <- mean(df$F0_Regression_Normalized) - 1*sd(df$F0_Regression_Normalized)

# Define a sequence of values for log_cum
log_cum_seq <- seq(from = min(df$log_cum), to = max(df$log_cum), length.out = 100)

# Create new data for both scenarios
new_data_low_F1 <- data.frame(
  F1_Regression_Normalized = constant_F1_low,
  F0_Regression_Normalized = mean(df$F0_Regression_Normalized, na.rm=T),
  log_cum = log_cum_seq,
  frt = "bck",
  CHILD = rep(unique(df$CHILD), each = length(log_cum_seq))
)

new_data_high_F0 <- data.frame(
  F0_Regression_Normalized = constant_F0_high,
  F1_Regression_Normalized = mean(df$F1_Regression_Normalized, na.rm=T),
  log_cum = log_cum_seq,
    frt = "bck",
  CHILD = rep(unique(df$CHILD), each = length(log_cum_seq))
)

new_data_high_F1 <- data.frame(
  F1_Regression_Normalized = constant_F1_high,
  F0_Regression_Normalized = mean(df$F0_Regression_Normalized, na.rm=T),
  log_cum = log_cum_seq,
    frt = "bck",
  CHILD = rep(unique(df$CHILD), each = length(log_cum_seq))
)

new_data_low_F0 <- data.frame(
  F0_Regression_Normalized = constant_F0_low,
  F1_Regression_Normalized = mean(df$F1_Regression_Normalized, na.rm=T),
  log_cum = log_cum_seq,
    frt = "bck",
  CHILD = rep(unique(df$CHILD), each = length(log_cum_seq))
)

# Predict probabilities for the new data
predicted_probs_low_F1 <- predict(final_model, newdata = new_data_low_F1, type = "response")
predicted_probs_high_F0 <- predict(final_model, newdata = new_data_high_F0, type = "response")
predicted_probs_high_F1 <- predict(final_model, newdata = new_data_high_F1, type = "response")
predicted_probs_low_F0 <- predict(final_model, newdata = new_data_low_F0, type = "response")

# Combine log_cum_seq, predicted_probs, and CHILD into data frames
predicted_df_low_F1 <- data.frame(
  log_cum = rep(log_cum_seq, length(unique(df$CHILD))),
  predicted_prob = predicted_probs_low_F1,
  CHILD = rep(unique(df$CHILD), each = length(log_cum_seq))
)

predicted_df_high_F0 <- data.frame(
  log_cum = rep(log_cum_seq, length(unique(df$CHILD))),
  predicted_prob = predicted_probs_high_F0,
  CHILD = rep(unique(df$CHILD), each = length(log_cum_seq))
)

predicted_df_high_F1 <- data.frame(
  log_cum = rep(log_cum_seq, length(unique(df$CHILD))),
  predicted_prob = predicted_probs_high_F1,
  CHILD = rep(unique(df$CHILD), each = length(log_cum_seq))
)

predicted_df_low_F0 <- data.frame(
  log_cum = rep(log_cum_seq, length(unique(df$CHILD))),
  predicted_prob = predicted_probs_low_F0,
  CHILD = rep(unique(df$CHILD), each = length(log_cum_seq))
)

# Calculate average predicted probabilities
average_prob_low_F1 <- predicted_df_low_F1 %>%
  group_by(log_cum) %>%
  summarize(average_prob = mean(predicted_prob))

average_prob_high_F0 <- predicted_df_high_F0 %>%
  group_by(log_cum) %>%
  summarize(average_prob = mean(predicted_prob))

average_prob_high_F1 <- predicted_df_high_F1 %>%
  group_by(log_cum) %>%
  summarize(average_prob = mean(predicted_prob))

average_prob_low_F0 <- predicted_df_low_F0 %>%
  group_by(log_cum) %>%
  summarize(average_prob = mean(predicted_prob))

# Plot the lines
plot1 <- ggplot(predicted_df_low_F1) +
  geom_line(aes(x = log_cum, y = predicted_prob, group = CHILD), color = "orange", alpha = 0.2) +
  geom_line(data = average_prob_low_F1, aes(x = log_cum, y = average_prob), color = "orange", size = 2) +
  labs(x = "", y = "", title = "a low F1 (mean-SD)") +
     scale_y_continuous(labels = percent_format(), limits = c(0, 1)) + 
  theme_minimal()

plot2 <- ggplot(predicted_df_high_F0) +
  geom_line(aes(x = log_cum, y = predicted_prob, group = CHILD), color = "purple", alpha = 0.2) +
  geom_line(data = average_prob_high_F0, aes(x = log_cum, y = average_prob), color = "purple", size = 2) +
  labs(x = "", y = "", title = "a high F0 (mean+SD)") +
     scale_y_continuous(labels = percent_format(), limits = c(0, 1)) + 
  theme_minimal()

plot3 <- ggplot(predicted_df_high_F1) +
  geom_line(aes(x = log_cum, y = predicted_prob, group = CHILD), color = "orange", alpha = 0.2) +
  geom_line(data = average_prob_high_F1, aes(x = log_cum, y = average_prob), color = "orange", size = 2) +
  labs(x = "log (cumulative vocabulary)", y = "", title = "a high F1 (mean+SD)") +
     scale_y_continuous(labels = percent_format(), limits = c(0, 1)) + 
  theme_minimal()

plot4 <- ggplot(predicted_df_low_F0) +
  geom_line(aes(x = log_cum, y = predicted_prob, group = CHILD), color = "purple", alpha = 0.2) +
  geom_line(data = average_prob_low_F0, aes(x = log_cum, y = average_prob), color = "purple", size = 2) +
  labs(x = "log (cumulative vocabulary)", y = "", title = "a low F0 (mean-SD)") +
     scale_y_continuous(labels = percent_format(), limits = c(0, 1)) + 
  theme_minimal()

# Combine plots
combined_plots <- grid.arrange(plot1, plot2, plot3, plot4, ncol = 2)

# Add overall title
overall_title <- ggdraw() + 
  draw_label("Individual trajectories for the effect of", size = 15, fontface = "bold", x = 0.5)

# Add single y-axis label
y_label <- ggdraw() + 
  draw_label("Predicted probability of a vowel being high", size = 12, angle = 90, 
             x = 0.015, y = -1, vjust = 0.5, hjust = 1)
# Print combined plot with overall title
grid.arrange(overall_title, y_label, combined_plots, ncol = 1, heights = c(0.05, 0.05, 0.9))

save.image(file = "analysis_F1_F0.RData")