---
title: "rMA Thesis analysis scripts"
author: "Patrik Hrabánek"
date: "2024-05-23"
output:
  html_document: default
  pdf_document: default
---

For my thesis, I have data from a perception experiment using an ABX task. The RQ-answering variable I was interested in is whether Korean listeners perceive an artificially created creaky lenis token as belonging to the 'fortis' category or to the 'lenis' category. My data structure is as shown in "table" below.


Import libraries:
```{r}
library(dplyr)
library(minqa)
options(width=200)
```

Import the full datasheet. The name of the .txt file corresponds to the "fortis" responses being coded as 1 and "lenis" responses being coded 0:
```{r}
table <- read.delim("anonymizedData.tsv", stringsAsFactors = TRUE)
head(table)
```
I have a lot of fillers, which do not tell me anything about the RQ, so I filter them out:

```{r}
exp_tokens <- filter(table, Type == "Exp")
head(exp_tokens)
```

The key variable here is the 'fxlChoice' column. A value of 'f' means the listener said that sound X (from the ABX set) is closer to the sound that was a 'fortis' sound (always sound A). Analogously, a value of 'l' means that the listener said that sound X was closer to the 'lenis' sound (always sound B). The column 'Creak' signifies if natural creakiness was present in the A sound. The column 'Word' contains the two words which we used in the experimental condition (tajut or tasha).

Each ABX set was played four times in two orders: ABX and BAX. This is signified in the 'Order' column. 

For the demographic questions, we have Age (continuous variable), Gender (M or F), Region (6 levels), SeoulYears (whether the participant spent 2 or more years in Seoul; Yes or No), and Mandarin (knowledge of Mandarin; Yes or No). The important thing about Region is: We expected the Seoul region to score differently than the other regions. BUT, the "Gyeonggi" region can still be considered to be a part of "The greater Seoul" and the dialect there is not known to differ very much from the Seoul dialect. Therefore, we decided to group Regions into three levels in the column SeoulRegion: Seoul, Gyeonggi, and Other (all other regions). Other columns are irrelevant at this point.

```{r}
#Gender
contrast <- cbind (c(-1/2, +1/2))
colnames (contrast) <- c("-F+M")
contrasts (exp_tokens$Gender) <- contrast
contrasts (exp_tokens$Gender)

#SeoulYears
contrast <- cbind (c(-1/2, +1/2))
colnames (contrast) <- c("-N+Y")
contrasts (exp_tokens$SeoulYears) <- contrast
contrasts (exp_tokens$SeoulYears)

#Mandarin
contrast <- cbind (c(-1/2, +1/2))
colnames (contrast) <- c("-N+Y")
contrasts (exp_tokens$Mandarin) <- contrast
contrasts (exp_tokens$Mandarin)

#Creak
contrast <- cbind (c(-1/2, +1/2))
colnames (contrast) <- c("-N+Y")
contrasts (exp_tokens$Creak) <- contrast
contrasts (exp_tokens$Creak)

#Order
contrast <- cbind (c(-1/2, +1/2))
colnames (contrast) <- c("-NORMAL+REVERSED")
contrasts (exp_tokens$Order) <- contrast
contrasts (exp_tokens$Order)

#Word
contrast <- cbind (c(-1/2, +1/2))
colnames (contrast) <- c("-taju+tasha")
contrasts (exp_tokens$Word) <- contrast
contrasts (exp_tokens$Word)

#SeoulRegion
levels(exp_tokens$SeoulRegion)
contrast <- cbind (c(0, +0.5, -0.5), c(+2/3, -1/3, -1/3))
colnames (contrast) <- c("-Seoul+Other", "-SeoulRest+Gyeonggi")
contrasts (exp_tokens$SeoulRegion) <- contrast
contrasts (exp_tokens$SeoulRegion)
```

I looked closer at my data and found that Age and Region are strongly correlated. Therefore, I chose to omit Age in the final model, as I was more interested in the regional differences than age differences.

```{r}

#Set age groups
age_breaks <- c(20, 30, 40, 50, 60)
age_labels <- c("21-30", "31-40", "41-50", "51-60")

# New column with age groups
exp_tokens$AgeGroup <- cut(exp_tokens$Age, 
                           breaks = age_breaks, 
                           labels = age_labels, 
                           right = FALSE)

head(exp_tokens)


#Contingency table
contingency_table <- table(exp_tokens$Region, exp_tokens$AgeGroup)
contingency_table

```

Setting "bobyqa" as the optimizer to help issues with model convergence (as per Titia's advice):
```{r}
library(lme4)
ctrl <- glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5))
```


```{r}
library (lme4)
model_glm <- lme4::glmer (fxlChoice ~ Creak * (Order + SeoulRegion + Mandarin + Gender + Word) + (Creak * (Order + Word) | Participant), control = ctrl, family=binomial, data=exp_tokens)
summary (model_glm)
model_glm <- lme4::glmer (fxlChoice ~ Creak * (Order + SeoulRegion + Mandarin + Gender + Word) + (Creak * Order | Participant), control = ctrl, family=binomial, data=exp_tokens)
summary (model_glm)
model_glm <- lme4::glmer (fxlChoice ~ Creak * (Order + SeoulRegion + Mandarin + Gender + Word) + (Order | Participant), control = ctrl, family=binomial, data=exp_tokens)
summary (model_glm)
model_glm <- lme4::glmer (fxlChoice ~ Creak * (Order + SeoulRegion + Mandarin + Gender + Word) + (Creak | Participant), control = ctrl, family=binomial, data=exp_tokens)
summary (model_glm)
model_glm <- lme4::glmer (fxlChoice ~ Creak * (Order + SeoulRegion + Mandarin + Gender + Word) + (Order + Word | Participant), control = ctrl, family=binomial, data=exp_tokens)
summary (model_glm)
model_glm <- lme4::glmer (fxlChoice ~ Creak * (Order + SeoulRegion + Mandarin + Gender + Word) + (Word | Participant), control = ctrl, family=binomial, data=exp_tokens)
summary (model_glm)
#final model
model_glm <- lme4::glmer (fxlChoice ~ Creak * (Order + SeoulRegion + Mandarin + Gender + Word) + (1 | Participant), control = ctrl, family=binomial, data=exp_tokens)
summary (model_glm)
```


Mean and SD calculation of fortis/lenis tokens:

VOT:
```{r}
fortis_tasha <- c(8.4, 12.5, 12.6, 22.4, 14.1, 14.2, 6.6, 6.4)
lenis_tasha <- c(75.8, 70.4, 74.7, 67.6, 60, 66)
fortis_tajut <- c(9, 9.7, 9.7, 10.8, 8.7, 10.6, 8, 8.2)
lenis_tajut <- c(63.4, 64, 64.2, 44.3, 58.9, 53, 64, 60)


mean(lenis_tasha)
sd(lenis_tasha)

mean(lenis_tajut)
sd(lenis_tajut)
```

F0:
```{r}
fortis_f0_tajut <- c(225, 235, 234, 218, 231, 233, 228, 237)
lenis_f0_tajut <- c(244, 248, 250, 236, 239, 249, 238, 245)

fortis_f0_tasha <- c(213, 219, 222, 215, 227, 229, 232, 232)
lenis_f0_tasha <- c(238, 220, 224, 230, 224, 220)

mean(lenis_f0_tasha)
sd(lenis_f0_tasha)


mean(lenis_f0_tajut)
sd(lenis_f0_tajut)
```

Transforming log-odds to odds + lower and upper 95% confidence intervals:
```{r}
odds <- exp(2.58)
upper_ci <- exp(1.98218 + 2 * 0.19321)
lower_ci <- exp(1.98218 - 2 * 0.19321)

odds
lower_ci
upper_ci
```