---
title: "checkedness_experiment"
author: "Chengjia Ye"
date: "2023-12-01"
output: html_document
---

# 1. Study introduction

The present study examines whether phonation type, vowel quality and duration play a role in perception of checked tones in Shanghai Chinese. A three-alternative forced-choice experiment was conducted in which 40 native listeners chose among three Chinese characters representing an open, a nasal-ending and a glottal-ending syllable according to the word they heard. Results showed main effects of phonation type, vowel quality and duration as well as the interaction between the former two with respect to both decisions and reaction time (RT). Two-way interactions between duration and the other two cues, however, could only be seen in RT. We therefore tentatively conclude that tenser phonation may be the primary perceptual cue to checked tones while more central vowel space and shorter duration are secondary, in line with findings in other tonal languages like Vietnamese and Cantonese. Additionally, we also found that the low vowel pair /ɐ/ and /ɑ/ and high register tones tended to bias perception towards checked tones. 

This .Rmd file documents the data processing, statistical analysis and plotting of this study.

# 2. Load packages

```{r load packages}
library(lme4)
library(lmerTest)
library(ggplot2)
library(forcats)
library(viridis)
library(hrbrthemes)
library(tidyr)
library(ggsignif)
library(emmeans)
options(width=200)
```

# 3. Plotting all responses

## 3.1 Load data

```{r load data}
table.dt <- read.delim("results.tsv",  stringsAsFactors=TRUE)
table.dt$group <- as.factor(table.dt$group)
table.dt <- table.dt[,-4:-5]  #remove district and trial number
```

## 3.2 Count responses

```{r count responses}
participant <- levels(table.dt$participant)
dt.nasalph <- table.dt[table.dt$phonation == 'nasal',]
dt.openph <- table.dt[table.dt$phonation == 'unchecked',]
dt.glottalph <- table.dt[table.dt$phonation == 'checked',]
cnt <- 1
nasal <- vector(length = 40 * 3)
open <- vector(length = 40 * 3)
glottal <- vector(length = 40 * 3)
for (part in participant) {
  dt.part <- dt.openph[dt.openph$participant == part,]
  nasal[cnt] <- nrow(dt.part[dt.part$decision == "s",])
  open[cnt] <- nrow(dt.part[dt.part$decision == "a",])
  glottal[cnt] <- nrow(dt.part[dt.part$decision == "d",])
  cnt <- cnt+1
}
for (part in participant) {
  dt.part <- dt.nasalph[dt.nasalph$participant == part,]
  nasal[cnt] <- nrow(dt.part[dt.part$decision == "s",])
  open[cnt] <- nrow(dt.part[dt.part$decision == "a",])
  glottal[cnt] <- nrow(dt.part[dt.part$decision == "d",])
  cnt <- cnt+1
}
for (part in participant) {
  dt.part <- dt.glottalph[dt.glottalph$participant == part,]
  nasal[cnt] <- nrow(dt.part[dt.part$decision == "s",])
  open[cnt] <- nrow(dt.part[dt.part$decision == "a",])
  glottal[cnt] <- nrow(dt.part[dt.part$decision == "d",])
  cnt <- cnt+1
}
```

## 3.3 Transfer count to percentage

```{r count to percentage}
# Transfer count to percentage
new_frame <- data.frame(participant = 1:40, open = c(open[1:40]/0.64, open[41:80]/0.32, open[81:120]/0.64), nasal = c(nasal[1:40]/0.64, nasal[41:80]/0.32, nasal[81:120]/0.64), glottal = c(glottal[1:40]/0.64, glottal[41:80]/0.32, glottal[81:120]/0.64), stimuli_type = rep(c("open", "nasal", "glottal"), each = 40))

new_frame_long <- gather(new_frame[,-1], key = "perceived_syllable", value = "percentage", -stimuli_type)
new_frame_long$stimuli_type <- factor(new_frame_long$stimuli_type, levels = c("open", "nasal", "glottal"))
new_frame_long$perceived_syllable <- factor(new_frame_long$perceived_syllable, levels = c("open", "nasal", "glottal"))
```

## 3.4 Plot

```{r plot all data}
fig_2 <- ggplot(new_frame_long, aes(x = stimuli_type, 
                           y= percentage, fill = perceived_syllable)) + ylab("Percentage of responses (%)") +
  scale_x_discrete(labels = c("[−tense] [−nasal]","[−tense] [+nasal]","[+tense] [−nasal]")) +
  scale_y_continuous(limits=c(0,100), breaks=seq(0,100,20), expand = c(0, 0)) +
  labs(fill = "Perceived syllable") + 
  geom_boxplot(outlier.size = 0.5, notch = FALSE, linewidth = 0.35, fatten = 2) +
  scale_fill_viridis(discrete = TRUE, alpha=0.5, option="D", labels = c("open", "nasal-ending", "glottal-ending")) +
  theme_ipsum() +
  theme(
      legend.position = "bottom",
      axis.title.x  = element_text(vjust = 0.5, hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.title.y  = element_text(hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 10.5),
      axis.text.y   = element_text(family = "serif", size = 10.5),
      legend.title  = element_text(family = "serif", face="bold",size = 12),
      legend.text = element_text(family = "serif", size = 10.5),
      legend.box.spacing = margin(0.25),
      plot.margin = unit(c(0.25,0.75,0.25,0.75), units = "cm")
    ) + xlab("Phonation type and final nasality as auditory cues in stimuli")
      #+ xlab("Stimuli sorted by two auditory cues: phonation and nasality")
ggsave("Figure_2.png", fig_2, width = 132, height = 72, units = "mm", dpi = 1000)
fig_2
```

# 4. Responses

## 4.1 Adjust data

```{r adjust data for response analysis}
#remove responses related to "nasal"
table.dt <- table.dt[table.dt$decision != "s" & table.dt$phonation != "nasal", ] 
table.dt <- droplevels(table.dt)
```

## 4.2 Contrast coding

```{r contrast coding}
# Duration: orthogonal polynomial coding
#contrasts(table.dt$duration) <- contr.poly(4)
contrast.duration <- cbind(c(-1/2,-1/6,1/6,1/2),c(0.5,-0.5,-0.5,0.5),c(-1/6,1/2,-1/2,1/6))
colnames(contrast.duration) <- c("L","Q","C")
contrasts(table.dt$duration) <- contrast.duration
contrasts(table.dt$duration)

# Vowel quality
contrast.vowel_quality <- cbind (c(+0.5,-0.5))
colnames(contrast.vowel_quality) <- c("+closed-open")
contrasts(table.dt$vowel_quality) <- contrast.vowel_quality
contrasts(table.dt$vowel_quality)

# Vowel pair
contrast.vowel_pair <- cbind (c(+0.5,-0.5))
colnames(contrast.vowel_pair) <- c("+a-o")
contrasts(table.dt$vowel_pair) <- contrast.vowel_pair
contrasts(table.dt$vowel_pair)

# Phonation type
contrast.phonation <- cbind (c(+0.5,-0.5))
colnames(contrast.phonation) <- c("+checked-unchecked")
contrasts(table.dt$phonation) <- contrast.phonation
contrasts(table.dt$phonation)

# Tone register
contrast.tone <- cbind (c(-0.5,+0.5))
colnames(contrast.tone) <- c("-high+low")
contrasts(table.dt$tone) <- contrast.tone
contrasts(table.dt$tone)
```

What about age? If we centre all our binary predictors, we should also centre Age.
Otherwise, the *p*-values for all predictors not involving Age will be for newborns.

```{r}
mean (table.dt$age)
sd (table.dt$age)
```

So let’ centre Age:
```{r}
table.dt$age_centred <- (table.dt$age - 27.0) / 10.0
```
The *p*-values for all predictors not involving age_centred will then be for 27-year-olds,
and the estimate of the effect of age_centred will be in units of "per decade".

## 4.3 Models of decisions (glmer)

```{r glmer model}
mdl.dcsion <- lme4::glmer(decision ~ vowel_pair + tone + vowel_quality*duration*phonation + (1|participant), data=table.dt, family=binomial)
summary(mdl.dcsion)
```
The model doesn't converge if Age is included:
```{r glmer model with age}
mdl.dcsion_age <- lme4::glmer(decision ~ age + vowel_pair + tone + vowel_quality*duration*phonation + (1|participant), data=table.dt, family=binomial)
summary(mdl.dcsion_age)
```
The model doesn't even converge if `age_centred` is included:
```{r glmer model with age_centred}
mdl.dcsion_age <- lme4::glmer(decision ~ age_centred + vowel_pair + tone + vowel_quality*duration*phonation + (1|participant), data=table.dt, family=binomial)
summary(mdl.dcsion_age)
```
Sometimes it helps if only `bobyqa` is used as the optimizer:
```{r}
ctrl <- glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5))
mdl.dcsion_age <- lme4::glmer(decision ~ age_centred + vowel_pair + tone + vowel_quality*duration*phonation + (1|participant), data=table.dt, family=binomial, control=ctrl)
summary(mdl.dcsion_age)
```
No problem, but the centring was crucial:
```{r}
mdl.dcsion_age <- lme4::glmer(decision ~ age + vowel_pair + tone + vowel_quality*duration*phonation + (1|participant), data=table.dt, family=binomial, control=ctrl)
summary(mdl.dcsion_age)
```

Thus, the age model converges only if we use `bobyqa` *and* centre the Age predictor!

Now, the *p*-value for age is around 0.70, so we can ignore this predictor from now on.

The maximal model would be:

```{r maximal model}
mdl.dcsion_maximal <- lme4::glmer(decision ~ vowel_pair*tone*vowel_quality*duration*phonation + (vowel_pair*tone*vowel_quality*duration*phonation|participant), data=table.dt, family=binomial, control=ctrl)
summary(mdl.dcsion_maximal)
```

## 4.4 Calculating odds ratio

```{r calculating odds ratio}
# odds calculating function
odds_calculating <- function(row, variable_name) {
  estimate <- abs(summary(mdl.dcsion)$coefficients[row,1])
  sd <- abs(summary(mdl.dcsion)$coefficients[row,2])
  exp_estimate <- round(exp(estimate),2)
  exp_interval <- c(round(exp(estimate - 1.96 * sd),2), round(exp(estimate + 1.96 * sd),2))
  z <- round(summary(mdl.dcsion)$coefficients[row,3],3)
  p <- summary(mdl.dcsion)$coefficients[row,4]
  print(c(variable_name, exp_estimate, exp_interval, z, p))
}

# vowel pair
odds_calculating(2, "vowel pair")

# tone
odds_calculating(3, "tone register")

# vowel quality
odds_calculating(4, "vowel quality")

# duration.L
odds_calculating(5, "duration linear")

#phonation type
odds_calculating(8, "phonation type")

# vowel quality and phonation type
odds_calculating(12, "vowel quality and phonation")
```

## 4.5 Plotting the interaction

```{r interaction plot}
# Load data and create vectors to accommodate data
cnt <- 1
open.ckd <- vector(length = 40)
closed.ckd <- vector(length = 40)
open.unckd <- vector(length = 40)
closed.unckd <- vector(length = 40)
for (part in participant) {
  dt.part <- table.dt[table.dt$participant == part,]
  open.ckd[cnt] <- nrow(dt.part[dt.part$phonation == 'checked' & dt.part$vowel_quality == "open" & dt.part$decision == "d",])
  closed.ckd[cnt] <- nrow(dt.part[dt.part$phonation == 'checked' & dt.part$vowel_quality == "closed"& dt.part$decision == "d",])
  open.unckd[cnt] <- nrow(dt.part[dt.part$phonation == 'unchecked' & dt.part$vowel_quality == "open"& dt.part$decision == "d",])
  closed.unckd[cnt] <- nrow(dt.part[dt.part$phonation == 'unchecked' & dt.part$vowel_quality == "closed"& dt.part$decision == "d",])
  cnt <- cnt+1
}

# gather data together
DF <- data.frame(responses = (c(c(open.ckd, closed.ckd),c(open.unckd,closed.unckd))/0.32), phonation = rep(c("tense", "modal"), each = 80), vowel_quality = rep(rep(c("unchecked","checked"),each = 40), 2), stringsAsFactors = FALSE)

DF$phonation <- factor(DF$phonation, levels = c("tense","modal"))
DF$responses <- DF$responses * 0.32

# plot
fig_3 <- ggplot(DF, aes(x = phonation,
                           y= responses, fill = vowel_quality)) + ylab("Stimuli perceived as checked tones") +
  scale_y_continuous(limits=c(0,35.5), breaks=seq(0,32,8), expand = c(0, 0)) +
  labs(fill = "Vowel quality") + 
  geom_boxplot(outlier.size = 0.75, notch = FALSE, linewidth = 0.35, fatten = 2) +
  #geom_violin() +
  scale_fill_viridis(discrete = TRUE, alpha=0.5, option="D", labels = c("checked\n/o/ and /ɐ/","unchecked\n/ɔ/ and /ɑ/")) +
  theme_ipsum() +
  theme(
      #legend.position = c(0.95,0.9),
      legend.position = "right",
      axis.title.x  = element_text(vjust = 0.5, hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.title.y  = element_text(hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 10.5),
      axis.text.y   = element_text(family = "serif", size = 10.5),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", size = 10.5, 
      margin = margin(b = 10)),
      plot.margin = unit(c(0.15,0.75,0.15,0.75), units = "cm")
    ) + xlab("Phonation type") +
   geom_signif(comparisons = list(c("tense", "modal")), map_signif_level = TRUE, textsize = 4, vjust = 0.1)
ggsave("Figure_3.png", fig_3, width = 130, height = 75, units = "mm", dpi = 1000)
fig_3
```

```{r plotting the fitted model}
predicted_values <- predict(mdl.dcsion, newdata = table.dt, type = "response")
#interaction.plot(table.dt$phonation, table.dt$vowel_quality, response = predicted_values, type = "b", col = c("red", "blue"))

ggplot(table.dt, aes(x = table.dt$phonation, y = predicted_values, color = factor( table.dt$vowel_quality))) + ylab("mean percentage of stimuli perceived as checked tones") +
  geom_point(position = position_dodge(width = 0.2)) +
  #geom_violin() +
  geom_line(aes(group = factor(table.dt$vowel_quality)), position = position_dodge(width = 0.2)) +
  labs(color = "Vowel quality") +
  scale_fill_viridis(discrete = TRUE, alpha=0.5, option="D") +
  theme_ipsum()+
  theme(
      #legend.position = c(0.95,0.9),
      legend.position = "right",
      axis.title.x  = element_text(vjust = 0.5, hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.title.y  = element_text(hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 10.5),
      axis.text.y   = element_text(family = "serif", size = 10.5),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", size = 10.5, margin = margin(b = 10))
    ) + xlab("phonation type") +
   geom_signif(comparisons = list(c("checked", "unchecked")), map_signif_level = TRUE, textsize = 4, vjust = 0.1)
```

```{r posthoc}
summary(pairs(emmeans(mdl.dcsion, ~ vowel_quality | phonation), adjust = "bonferroni"))
```

# 5. Reaction time

## 5.1 Data preparation

```{r RT}
# reload data
table.dt2 <- read.delim("results.tsv",  stringsAsFactors=TRUE)
table.dt2 <- table.dt2[,-4:-5]
table.dt2$group <- as.factor(table.dt2$group)

# remove responses as nasal
table.dt2 <- table.dt2[table.dt2$phonation != "nasal", ]
(nrow(table.dt2))
table.dt2 <- droplevels(table.dt2)

#remove responses with RT longer than 5125 ms
table.dt2 <- table.dt2[table.dt2$rt <= 5125, ]
(nrow(table.dt2))

# contrast coding
#contrasts(table.dt2$duration) <- contr.poly(4)
#contrasts(table.dt2$duration)
contrasts(table.dt2$duration) <- contrast.duration
contrasts(table.dt2$duration)

contrasts(table.dt2$vowel_quality) <- contrast.vowel_quality
contrasts(table.dt2$vowel_pair) <- contrast.vowel_pair
contrasts(table.dt2$phonation) <- contrast.phonation
contrasts(table.dt2$tone) <- contrast.tone
```

## 5.2 Statistical model

```{r model RT}
mdl.rt <- lmerTest::lmer(rt ~ vowel_pair + tone + vowel_quality*duration*phonation + (1|participant), data=table.dt2)
summary(mdl.rt)
```
```{r model RT with age}
mdl.rt_age <- lmerTest::lmer(rt ~ age + vowel_pair + tone + vowel_quality*duration*phonation + (1|participant), data=table.dt2)
summary(mdl.rt_age)
```

No effect of age on reaction time either.

```{r plotting the interaction}
ggplot(table.dt2, aes(x = phonation,
                           y= rt, fill = vowel_quality)) + 
  geom_boxplot(na.rm = TRUE) + 
  scale_y_continuous(limits=c(500,2500), breaks=seq(500,2500,500), expand = c(0, 0)) +
  ylab("Reaction time (ms)") +
  scale_fill_viridis(discrete = TRUE, alpha=0.5, option="D") +
  labs(fill = "Vowel quality") + 
  theme_ipsum() +
  theme(
      legend.position = "right",
      axis.title.x  = element_text(vjust = 0.5, hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.title.y  = element_text(hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 10.5),
      axis.text.y   = element_text(family = "serif", size = 10.5),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", size = 10.5, margin = margin(b = 10))
    )
```

```{r plotting the fitted model 2}
predicted_values <- predict(mdl.rt, newdata = table.dt2, type = "response")
ggplot(table.dt2, aes(x = table.dt2$phonation, y = predicted_values, color = factor( table.dt2$vowel_quality))) + ylab("Reaction time (ms)") +
  geom_point(position = position_dodge(width = 0.2)) +
  #geom_violin() +
  geom_line(aes(group = factor(table.dt2$vowel_quality)), position = position_dodge(width = 0.2)) +
  labs(color = "Vowel quality") +
  scale_fill_viridis(discrete = TRUE, alpha=0.5, option="D") +
  theme_ipsum()+
  theme(
      #legend.position = c(0.95,0.9),
      legend.position = "right",
      axis.title.x  = element_text(vjust = 0.5, hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.title.y  = element_text(hjust = 0.5, family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 10.5),
      axis.text.y   = element_text(family = "serif", size = 10.5),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", size = 10.5, margin = margin(b = 10))
    ) + xlab("phonation type") +
   geom_signif(comparisons = list(c("checked", "unchecked")), map_signif_level = TRUE, textsize = 4, vjust = 0.1)
```

## 5.3 Calculating intervals

```{r calculating intervals}
# function
interval_calculating <- function(row, variable_name) {
  estimate <- summary(mdl.rt)$coefficients[row,1]
  sd <- summary(mdl.rt)$coefficients[row,2]
  rounded_estimate <- round(estimate)
  interval <- c(round((estimate - 1.96 * sd),0), round((estimate + 1.96 * sd),0))
  df <- summary(mdl.rt)$coefficients[row,3]
  df <- round(df,2)
  t <- round(summary(mdl.rt)$coefficients[row,4],3)
  p <- summary(mdl.rt)$coefficients[row,5]
  return(c(variable_name, rounded_estimate, interval, df, t, p))
}

# vowel pair
interval_calculating(2, "vowel pair")

# tone
interval_calculating(3, "tone register")

# vowel quality
interval_calculating(4, "vowel quality")

# duration.L
interval_calculating(5, "duration linear")

#phonation type
interval_calculating(8, "phonation type")

# vowel quality and phonation type
interval_calculating(12, "vowel quality and phonation")

# vowel quality and duration.L
interval_calculating(9, "vowel quality and duration.L")

# phonation type and duration.L
interval_calculating(13, "phonation type and duration.L")
```
