# ------------------------------------------------------------------------------------------------------------------------
# This is a script to reproduce the findings and figures from
# 
# Kleiner, C., & Birkholz, P. (2023). Intrinsic velocity differences between larynx raising and larynx lowering. PloS One. 
#
# Input file: S1_File.csv
#
# Christian Kleiner, 25.01.2023
# ------------------------------------------------------------------------------------------------------------------------

# Load packages.
library(dplyr)      ## for data manipulation (using pipes with %>%)
library(ggplot2)    ## for graphic functions
library(lmerTest)   ## to fit linear mixed models
library(emmeans)    ## to better understand fitted statistical regression models: predictions & comparisons
library(ggpubr)     ## to combine multiple plots
library(report)     ## to generate standard report automatically for guidance
library(MuMIn)      ## to calculate conditional and marginal R2

# Load data.
data = read.csv("S1_File.csv")

# ------------------------------------------------------------------------------------------------------------------------
# Fit mixed-effects model using effect coding. 

# Aggregate data.
data %>%
  bind_rows(data, .id = "id") %>%
  mutate(v_cm_per_s = ifelse(id == 1, vUp_cm_per_s, vDown_cm_per_s)) %>%
  mutate(dir = factor(ifelse(id == 1, -0.5, 0.5))) %>%				  # effect-coded movement direction (up: -0.5, down: 0.5)
  mutate(speakingRate = factor(ifelse(speakingRate == "slow", -0.5, 0.5))) %>%  # effect-coded speaking rate (slow: -0.5, fast: 0.5)
  mutate(sex = factor(substr(subject, start=1, stop=1)=="F", -0.5, 0.5)) %>%    # effect-coded sex (female: -0.5, male: 0.5) 
  select(subject, sex, vowelPair, vowelPairOrder, speakingRate, measurementNumber, v_cm_per_s, dir) %>%
  group_by(subject, sex, vowelPair, speakingRate, dir) %>%
  summarise(v_cm_per_s = median(v_cm_per_s, na.rm = TRUE)) %>%  
  arrange(subject, speakingRate, vowelPair, dir) -> data_temp

# Fit model.
lmm <- lmer(log(v_cm_per_s) ~ dir + sex/dir + speakingRate/dir + (speakingRate+dir|subject) + (1|vowelPair), data_temp)
summary(lmm)

# Diagnostic plots to validate homoscedasticity and normality of residuals, and normality of random effect
par(mfrow = c(1,3))
resid <- resid(lmm, type = "pearson")
ranef <- ranef(lmm)[["subject"]][[1]]
plot(resid ~ fitted(lmm), main = "Residual vs fitted"); abline(h = 0, lty = "dashed", col = "red")
qqnorm(y = resid, main = "Residual QQ-plot"); qqline(y = resid)
qqnorm(y = ranef, main = "Random Effect QQ-plot"); qqline(y = ranef)

# Predicted velocity ratios (larynx lowering vs. larynx raising) depending on sex, with p values.
lmm %>% 
  emmeans(~dir|sex) %>% 
  contrast(method="revpairwise") %>% 
  test(null = 0, side = ">", type = "response", by = NULL)

# ------------------------------------------------------------------------------------------------------------------------
# Fig 3 from paper.

# Aggregate data.
data %>%
  bind_rows(data, .id = "id") %>%
  mutate(dir = ifelse(id == 1, "up", "down")) %>%
  mutate(speakingRate = factor(speakingRate, levels=c("slow", "fast"))) %>%
  mutate(sex = factor(substr(subject, start=1, stop=1), levels=c("F", "M"))) %>%
  mutate(dir = factor(dir, levels=c("up", "down"))) %>%
  mutate(v_cm_per_s = ifelse(id == 1, vUp_cm_per_s, vDown_cm_per_s)) %>%
  select(subject, sex, vowelPair, vowelPairOrder, speakingRate, measurementNumber, v_cm_per_s, dir) %>%
  group_by(subject, sex, vowelPair, speakingRate, dir) %>%
  summarise(v_cm_per_s = median(v_cm_per_s, na.rm = TRUE)) %>%  
  arrange(subject, speakingRate, vowelPair, dir) -> data_temp
data_temp %>%
  mutate(subject = ifelse(sex == "M", "M01-M10", "F01-F10")) %>%
  bind_rows(data_temp) -> data_temp

# Generate figure.
plt1 <- ggplot2::ggplot(data = data_temp, mapping = aes(x = subject, y = v_cm_per_s, fill = interaction(dir, speakingRate))) +
  geom_boxplot(lwd=0.1, outlier.size = 0.1) +
  labs(x = "Subject", y = "", fill = "") +
  scale_y_continuous(trans='log10') + 
  scale_fill_grey(start=1, end=0.7, labels=c('Upward at slow speaking rate', 'Downward at slow speaking rate', 'Upward at fast speaking rate', 'Downward at fast speaking rate')) + 
  theme_minimal() + 
  theme(
    legend.position = c(0.8, 0.2), 
    text=element_text(size=10, family="serif"), 
    axis.text.x = element_text(size = 8), 
    panel.grid = element_line(size = 0.1), 
    axis.ticks = element_line(colour = "black", size = 0.1),
    axis.text.y = element_text(size=8),
    panel.border = element_rect(colour = "black", fill=NA, size=0.1))
plt2 <- ggplot2::ggplot(data = data_temp, mapping = aes(x = vowelPair, y = v_cm_per_s, fill = interaction(dir, speakingRate))) +
  geom_boxplot(lwd=0.1, outlier.size = 0.1) +
  labs(x = "Vowel pair", y = "Vertical larynx velocity [cm/s]", fill = "") + 
  scale_y_continuous(trans='log10') + 
  scale_fill_grey(start=1, end=0.7) + 
  theme_minimal() + 
  theme(
    legend.position = "none", 
    text=element_text(size=10, family="serif"), 
    axis.text.x = element_text(size = 8), 
    panel.grid = element_line(size = 0.1), 
    axis.ticks = element_line(colour = "black", size = 0.1),
    axis.text.y = element_text(size=8),
    panel.border = element_rect(colour = "black", fill=NA, size=0.1))
ggarrange(plt1, plt2, nrow=2, ncol=1)
ggsave("fig3.pdf", width = 19.05, height = 10, units = "cm")

# ------------------------------------------------------------------------------------------------------------------------
# Fig 4 from paper.

# Aggregate data.
data %>% 
  bind_rows(data, .id = "id") %>%
  mutate(dir = ifelse(id == 1, "up", "down")) %>%
  mutate(speakingRate = factor(speakingRate, levels=c("slow", "fast"))) %>%
  mutate(sex = factor(substr(subject, start=1, stop=1), levels=c("F", "M"))) %>%
  mutate(dir = factor(dir, levels=c("up", "down"))) %>%
  mutate(v_cm_per_s = ifelse(id == 1, vUp_cm_per_s, vDown_cm_per_s)) %>%
  group_by(subject, sex, vowelPair, vowelPairOrder, speakingRate) %>% 
  summarise(v_cm_per_s = median(v_cm_per_s, na.rm = TRUE), movementRange_cm = median(movementRange_cm, na.rm = TRUE)) -> data_temp

# Generate figure.
ggplot(data_temp, aes(x=movementRange_cm, y=v_cm_per_s, color = sex, shape=sex)) +
  geom_point(aes(size=sex), shape =16) +  
  scale_size_discrete(range = c(1,0.5)) + 
  scale_colour_grey(start=0.8, end=0) +
  labs(x = "", y = "", fill = "") + 
  theme_minimal() + 
  geom_smooth(method=lm, se = FALSE) + 
  theme(
    legend.position = "none",
    text=element_text(size=10, family="serif"), 
    axis.text.x = element_text(size = 8), 
    panel.grid = element_line(size = 0.1), 
    axis.ticks = element_line(colour = "black", size = 0.1),
    axis.text.y = element_text(size=8),
    panel.border = element_rect(colour = "black", fill=NA, size=0.1))
ggsave("fig4.pdf", width = 13, height = 5, units = "cm")

# ------------------------------------------------------------------------------------------------------------------------
# Fig 5 from paper.

# Aggregate data.
data %>% 
  bind_rows(data, .id = "id") %>%
  mutate(dir = ifelse(id == 1, "up", "down")) %>%
  mutate(speakingRate = factor(speakingRate, levels=c("slow", "fast"))) %>%
  mutate(sex = factor(substr(subject, start=1, stop=1), levels=c("F", "M"))) %>%
  mutate(dir = factor(dir, levels=c("up", "down"))) %>%
  mutate(v_cm_per_s = ifelse(id == 1, vUp_cm_per_s, vDown_cm_per_s)) -> data_temp
data_temp %>% 
  group_by(subject) %>% 
  summarise(v_cm_per_s = median(v_cm_per_s, na.rm = TRUE)) -> data_subject
data_temp%>% 
  group_by(subject) %>% 
  summarise(movementRange_cm = mean(movementRange_cm, na.rm = TRUE)) %>% 
  merge(data_subject, by="subject") -> data_subject
data_temp %>% group_by(subject, dir) %>% 
  summarise(v_cm_per_s = median(v_cm_per_s, na.rm = TRUE)) %>%
  summarise(v_down_up = last(v_cm_per_s) / first(v_cm_per_s)) %>%
  merge(data_subject, by="subject")-> data_subject
data_temp %>% 
  group_by(vowelPair) %>% 
  summarise(v_cm_per_s = median(v_cm_per_s, na.rm = TRUE)) -> data_vowelPair
data_temp %>% 
  group_by(vowelPair) %>% 
  summarise(movementRange_cm = mean(movementRange_cm, na.rm = TRUE)) %>% 
  merge(data_vowelPair, by="vowelPair") -> data_vowelPair

# Generate figure.
ggplot(data_subject, aes(x=movementRange_cm, y=v_cm_per_s)) +
  geom_point(color = "black", size = 2) +
  geom_text(aes(label=subject), size=2.5, family = "serif") +
  geom_smooth(se=FALSE, method="lm", color="black") + 
  geom_point(data = data_vowelPair, color = "#666666", shape=15, size = 2) +
  geom_text(data = data_vowelPair, aes(label=vowelPair), size=5, family = "serif") +
  geom_smooth(data = data_vowelPair, se=FALSE, method="lm", color="#666666") + 
  labs(x = "", y = "", fill = "") + 
  theme_minimal() + 
  theme(
    text=element_text(size=10, family="serif"), 
    axis.text.x = element_text(size = 8), 
    panel.grid = element_line(size = 0.1), 
    axis.ticks = element_line(colour = "black", size = 0.1),
    axis.text.y = element_text(size=8),
    panel.border = element_rect(colour = "black", fill=NA, size=0.1))
ggsave("fig5.pdf", width = 13, height = 5, units = "cm")

# ------------------------------------------------------------------------------------------------------------------------
# Additional anlyses.

# Significance and coefficients of linear relationships between mean vertical movement range and mean vertical velocity of 
# the larynx across the subjects and across the vowel pairs.
cor.test(data_subject$v_cm_per_s, data_subject$movementRange_cm)
cor.test(data_vowelPair$v_cm_per_s, data_vowelPair$movementRange_cm)
model <- lm(v_cm_per_s ~ movementRange_cm, data = data_subject)
summary(model)
model <- lm(v_cm_per_s ~ movementRange_cm, data = data_vowelPair)
summary(model)

# No significant difference in movement range between males and females.
t.test(data_subject$movementRange_cm[1:10], data_subject$movementRange_cm[11:20])

# Significance and coefficients of sex-dependent linear relationships between mean vertical movement range and vertical 
# velocity of the larynx.
data_female  <- data_temp[1:640,]
data_male <- data_temp[641:1280,]
cor.test(data_female$movementRange_cm, data_female$v_cm_per_s)
cor.test(data_male$movementRange_cm, data_male$v_cm_per_s)
model <- lm(v_cm_per_s ~ movementRange_cm, data = data_female)
summary(model)
model <- lm(v_cm_per_s ~ movementRange_cm, data = data_male)
summary(model)

# Non-significant sex-dependent difference in by-subject slopes (mean vertical movement range vs vertical velocity of larynx). 
# Linear slope is associated to spring-mass system with linear mass-normalized stiffness (see paper).
data_temp %>% 
  group_by(subject,sex) %>% 
  do(model = lm(v_cm_per_s ~ movementRange_cm, data = .)) %>% 
  mutate(slope = model$coef[2]) %>% 
  filter(subject!="M01") -> coef
# Significant sex-dependent slope
ggplot(coef, aes(x=sex, y=slope)) +
  geom_boxplot(outlier.shape = NA)
t.test(coef$slope[1:10], coef$slope[11:19])  

# Significance of linear by-subject slopes.
data_temp %>% 
  group_by(subject, sex) %>% 
  summarize(cor=cor(movementRange_cm, v_cm_per_s, use = "complete.obs"), p=cor.test(movementRange_cm, v_cm_per_s, use = "complete.obs")$p.value)

# No very high correlation between by-subject slopes and velocity ratios.
cor.test(coef$slope[1:10], data_subject$v_down_up[1:10])
cor.test(coef$slope[11:20], data_subject$v_down_up[11:20])