---
title: "QJEP_analysis"
author: "Leendert Plug"
date: "14 August 2023"
output:
  html_document:
    code_folding: hide
    df_print: paged
    theme: default
    toc: yes
  pdf_document:
    toc: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, include=FALSE}
#Package list
library(lme4)
library(lmerTest)
library(languageR)
library(partykit)
library(ggplot2)
library(diptest)
library(plyr)
library(pastecs)
library(cluster)
library(tidyverse)
library(factoextra)
library(randomForest)
library(ggeffects)
library(ordinal)
library(sjPlot)
library(devtools)
library(rstatix)
library(reshape2)
library(irr)
library(psych)
library(ggpubr)
library(car)
library(ggcorrplot)
library(glmmTMB)
library(emmeans)

#General plotting theme
theme_set(theme_bw(base_size = 24))
```
***

#Experiment 1

##Data structure

All of the results data are in *QJEP_Exp1_results.csv*. This includes responses to practice trials (block 'p'), which we can exclude for analysis purposes:

```{r}
data <- read.csv("QJEP_Exp1_results.csv",T)
datam <- data[data$block != "p",]
```

Note that Robert's original results spreadsheet had 'xx' for n.a. cells. I have made these empty cells. Robert also recorded survey results on the 'I saw' sentences which we changed to 'I see' as numbers prefaced by 'xx'. I have removed the 'xx', so these survey results (which were not exactly on the sentences we ended up with) are not flagged in this dataset.

Checking variable classes and fixing where relevant:

```{r}
sapply(datam, class)
datam$participant <- as.factor(datam$participant)
datam$rating <- as.numeric(as.character(datam$rating))
datam$trial <- as.numeric(as.character(datam$trial))
datam$block_trial <- as.numeric(as.character(datam$block_trial))
datam$Block <- as.factor(datam$block)
sapply(datam, class)
```

The crucial variable is *rating*, which runs between -500 and +500. Adding 500 to get rid of negatives. 

```{r}
datam$ratingm <- (datam$rating) + 500
min(datam$ratingm)
max(datam$ratingm)

datam$ratingp <- (datam$ratingm)/1000
min(datam$ratingp)
max(datam$ratingp)

datam$ratingp[datam$ratingp == 0] <- 0.0001
datam$ratingp[datam$ratingp == 1] <- 0.9999
min(datam$ratingp)
max(datam$ratingp)
```

Also getting rid of negatives in *x_manipulation_diff* and *x_semantic_diff* (the difference measures for manipulation goodness and semantic goodness) -- making the scale centered around 100 for *x_manipulation_diff* and around 10 for *x_semantic_diff*. 

```{r}
datam$x_manipulation_diffm <- datam$x_manipulation_diff + 100
min(datam$x_manipulation_diffm)
max(datam$x_manipulation_diffm)

datam$x_semantic_diffm <- datam$x_semantic_diff + 10
min(datam$x_semantic_diffm)
max(datam$x_semantic_diffm)
```

Now separating filler and experimental trials:

```{r}
datame <- datam[datam$ortho_version != "filler",]
datamf <- datam[datam$ortho_version == "filler",]
```

##Descriptives

Ratings (on 0-1000 scale) do not appear normally distributed across filenames and participants: there is a tendency for participants to hear x as similar to b (the faster comparison sentence). Given this shape, a log transformation does not increase normality.  

```{r}
datame %>%
  ggplot(aes(x=ratingm)) + 
geom_histogram(color="black", lwd=1, binwidth = 125) + 
  xlab("Rating") + ylab("Count") +
  scale_x_continuous(breaks = seq(0,1000,250)) + expand_limits(x = 0) 
```

Trying square root transformation -- still .

```{r}
datame$ratingm1 <- sqrt(max(datame$ratingm + 1) - datame$ratingm)
datame %>%
  ggplot(aes(x=ratingm1)) + 
  geom_density(color="black", fill="grey") +
   xlab("Rating") 
```

##Modelling

We can proceed with the base model plus *block* and run the addition of goodness and orthography to this model. First adding *x_semantic_goodness* and *x_manip_goodness* -- this suggests as above that *x_semantic_goodness* should be kept. 


```{r}
mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block, data = datame) 
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + x_semantic_goodness, data = datame) 
anova(mod1, mod2)

mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + x_semantic_goodness, data = datame) 
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + x_semantic_goodness + x_manip_goodness, data = datame) 
anova(mod1, mod2)
```

On recommendation of a Speech Prosody reviewer I am also running these models as beta family models -- as the rating scale is bounded on both ends (as I understand it). To run these models using the glmmTMB package I have had to manually replace 0 and 1 values by 0.0001 and 0.9999 respectively; I don't know why but including the extremes generates error messages. At this step, these models show the same as the lmer ones. 

```{r}
mod1 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block, data = datame, family=beta_family(link="logit")) 
mod2 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + x_semantic_goodness, data = datame, family=beta_family(link="logit")) 
anova(mod1, mod2)

mod1 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + x_semantic_goodness, data = datame, family=beta_family(link="logit")) 
mod2 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + x_semantic_goodness + x_manip_goodness, data = datame, family=beta_family(link="logit")) 
anova(mod1, mod2)
```

Now adding *ortho_version* to expanded base model (with *block* and *x_semantic_goodness*). As above, this suggests that *ortho_version* is the stronger predictor and *x_semantic_goodness* adds no prediction to it.

```{r}
mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + x_semantic_goodness, data = datame) 
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + x_semantic_goodness + ortho_version, data = datame) 
anova(mod1, mod2)
summary(mod2)
```

Again the beta family model comparison shows the same. 

```{r}
mod1 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + x_semantic_goodness, data = datame, family=beta_family(link="logit")) 
mod2 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + x_semantic_goodness + ortho_version, data = datame, family=beta_family(link="logit")) 
anova(mod1, mod2)
summary(mod2)
```

Now testing for possibility of interaction between *ortho_version* and *block*. This significantly improves fit.

```{r}
mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version, data = datame) 
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
anova(mod1, mod2)
summary(mod2)
```

```{r}
mod1 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + ortho_version, data = datame, family=beta_family(link="logit")) 
mod2 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame, family=beta_family(link="logit")) 
anova(mod1, mod2)
summary(mod2)
```

Recall that participants were divided into two groups. Testing whether group assignment modulates the effect of *ortho_version* by adding an interaction term; no difference.

```{r}
mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version + group:ortho_version, data = datame)
anova(mod1, mod2)
```

```{r}
mod1 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame, family=beta_family(link="logit")) 
mod2 <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version + group:ortho_version, data = datame, family=beta_family(link="logit")) 
anova(mod1, mod2)
```

Similarly testing whether the *ortho_version* effect is modulated by the difference in semantic or manipulation goodness between the schwa and no-schwa interpretations. No evidence of this.

```{r}
mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version + x_semantic_diffm:ortho_version, data = datame)
anova(mod1, mod2)

mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version + x_manipulation_diffm:ortho_version, data = datame)
anova(mod1, mod2)
```

So the best model seems the following. In the resulting model, the intercept is the estimated response (close to zero on the original scale) in block 1 when the first presented is b and orthography is schwa-less. In block 2, and when the first presented is a, and when the orthography is schwaful, the rating is higher. The orthography effect varies significantly with block such that it is smaller in block 2. 

```{r}
mod <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
summary(mod)

#equivalent model, more efficient syntax

mod <- lmer(ratingm ~ (1|item) + (1|participant) + block * ortho_version, data = datame)
summary(mod)

emm <- emmeans(mod, pairwise ~ ortho_version | block, pbkrtest.limit = 5084)
emm

#glmmTMB model

modg <- glmmTMB(ratingp ~ (1|item) + (1|participant) + block * ortho_version, data = datame, family=beta_family(link="logit")) 
summary(modg)

```

Calculating estimated means and plotting some of them.

```{r}
mod <- lmer(ratingm ~ (1|item) + (1|participant) + Block * Orthography, data = datame)
summary(mod)

emm_options(pbkrtest.limit = 3640)
emmip(mod, Orthography ~ Block, CIs = TRUE, col = "black",
      linearg = list(linetype = "solid", lwd = 1), dotarg = list(size = 4.5), CIarg = list(lwd = 0.5, alpha = 1)) + 
  ggplot2::theme_bw(base_size = 20) +
    ggplot2::xlab("Block")+
    ggplot2::ylab("Estimated marginal mean")
```

Finally we can check we can check whether any other control variables have an impact. We can start with the final model in the previous section, i.e. with the design variables in to the extent that they elucidate the effect of *ortho_version*. None of the acoustic variables (*x_dur*, *x_intensity*, *F0.mean.hz*, *F0.kurt.hz*) emerges as significant. 

```{r}
mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version + log(x_dur), data = datame)
anova(mod1, mod2)

mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version + scale(x_intensity), data = datame)
anova(mod1, mod2)

mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version + scale(x_F0.mean.hz), data = datame)
anova(mod1, mod2)

mod1 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datame)
mod2 <- lmer(ratingm ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version + scale(x_F0.kurt.hz), data = datame)
anova(mod1, mod2)
```

Also checking that the optimal model for the full stimulus set looks pretty much the same for the smaller stimulus set that we used in Exp2. This is coded in variable *included_2b*.

```{r}
datasel <- datame[datame$included_2b == "yes",]

mod <- lmer(ratingp ~ (1|item) + (1|participant) + block + ortho_version + block:ortho_version, data = datasel)
summary(mod)
```


#Experiment 2

##Data structure

All of the results data are in *QJEP_Exp2_results.csv*. This does not include responses to practice trials. It does include responses to both filler and experimental trials. We separate these out here to create two data frames. We also check and fix variable classes.

```{r}
datab <- read.csv("QJEP_Exp2_results.csv",T)

sapply(datab, class)
datab$group <- as.factor(datab$group)
datab$trial <- as.numeric(as.character(datab$trial))
datab$stimulus <- as.factor(datab$stimulus)
datab$syll <- as.numeric(datab$syll)
datab$response2 <- as.factor(datab$response2)
datab$Block <- as.factor(datab$Block)
datab$Pair <- as.factor(datab$Pair)
sapply(datab, class)

datax <- datab[datab$type == "exp",]
dataf <- datab[datab$type == "fill",]
```

##Descriptives

The overall response distribution shows that 'same' responses are most common, and 'faster' responses are around twice as common as 'slower' responses. In other words, listeners often heard no difference between the pair members, and if they heard a difference they were biased towards hearing the second pair member as faster. This fits with the results of our previous experiments with a pairwise-comparison design, and others as well (and may be explainable in terms of a recency effect).

```{r}
datax$response1 <- factor(datax$response1, levels = c("slower", "same", "faster"))

ggplot(datax, aes(x=response1)) +
  geom_bar()
```

The *response1* variable codes how the second pair member is judged in comparison with the first. As the ambiguous forms occur in both first and second position, it is useful to recode so that the response records how the ambiguous form compares to its comparison sentence. This is *response2*, with '1' meaning slower, '2' same and '3' faster. To make the graphs easier to interpret, I've derived *response3* with 'a_slow', 'a_same' and 'a_fast' as the levels ('a' standing for 'ambiguous'). *response4* has 'slower', 'same' and 'faster' in the same configuration. 

As might be expected, 'a_same' is the majority response and 'a_fast' and 'a_slow' have quite similar counts with '3' predominating slightly. This is presumably because in about half the 'faster' and 'slower' cases above, the ambiguous sentence is in first position -- so about half of the 'faster' judgements for *response1* become 'a_slow' for *response3*, and about half of the (less common) 'slower' judgements for *response1* become 'a_fast' for *response3*. This redresses the balance between 'faster' and 'slower'.

```{r}
ggplot(datax, aes(x=response3)) +
  geom_bar()

datax$Response <- datax$response4
datax$Response <- factor(datax$Response, levels = c("slower", "same", "faster"))

ggplot(datax, aes(x=Response)) +
  geom_bar()
```

As indicated, the position of the ambiguous sentence varied. Plotting by *Pair_type* illustrates the overall bias of listeners towards hearing the second pair member as faster: when the ambiguous form is in first position, it is more often heard as slower than its comparison sentence; when it is in second position it is more often heard as faster.

```{r}
ggplot(datax, aes(x=Response)) +
  geom_bar(color="black", lwd=1) + 
  ylab("Count") +
  facet_wrap(~Pair)
```

##Modelling

We can first try the same modelling as on Exp1, treating *response* as a numerical variable (slower = 1, same = 2, faster = 3)

```{r}
datax$response2 <- as.numeric(datax$response2)
datax$trial <- as.numeric(datax$trial)
sapply(datax, class)
```

Base model:

```{r}
mod1 <- lmer(response2 ~ (1|item) + (1|participant), data = datax) 
summary(mod1)
```

Adding design features suggests that with *position* and *block* in the model, *trial* no longer adds significant prediction. 

```{r}
mod1 <- lmer(response2 ~ (1|item) + (1|participant), data = datax) 
mod2 <- lmer(response2 ~ (1|item) + (1|participant) + Pair, data = datax) 
anova(mod1, mod2)

mod1 <- lmer(response2 ~ (1|item) + (1|participant) + Pair, data = datax) 
mod2 <- lmer(response2 ~ (1|item) + (1|participant) + Pair + Block, data = datax) 
anova(mod1, mod2)

mod1 <- lmer(response2 ~ (1|item) + (1|participant) + Pair + Block, data = datax) 
mod2 <- lmer(response2 ~ (1|item) + (1|participant) + Pair + Block + trial, data = datax) 
anova(mod1, mod2)
```

Adding *orthography* to model with *position* and *block*: not significant as a main effect.

```{r}
mod1 <- lmer(response2 ~ (1|item) + (1|participant) + Pair + Block, data = datax) 
mod2 <- lmer(response2 ~ (1|item) + (1|participant) + Pair + Block + Orthography, data = datax) 
anova(mod1, mod2)
```

Adding *orthography* to model with *position* and *block* as an interaction with *position*: approaching significance if compared with the model with only *position* and *block* (two degrees of freedom, as both main effect and interaction effect get added), but adding interaction to main effects model only does yield a significant improvement of fit. Model summary gives significance for 'positiony:orthographyschwa' with a positive coefficient, which I take to reflect the descriptive observation that when the ambiguous form appears in second sentence, the 'schwa' condition yields a greater proportion of 'a_fast' responses than the 'noschwa' condition.

```{r}
mod1 <- lmer(response2 ~ (1|item) + (1|participant) + Pair + Block + Orthography, data = datax) 
mod2 <- lmer(response2 ~ (1|item) + (1|participant) + Block + Pair*Orthography, data = datax)
mod3 <- lmer(response2 ~ (1|item) + (1|participant) + Block + Pair + Orthography + Pair:Orthography, data = datax)
anova(mod1, mod2)
summary(mod2)

mod2 <- lmer(response2 ~ (1|item) + (1|participant) + Block + Pair*Orthography, data = datax)
mod4 <- lmer(response2 ~ (1|item) + (1|participant) + Block*Orthography + Pair*Orthography, data = datax)
anova(mod2, mod4)
```

Post-hoc comparisons 

```{r}
mod1 <- lmer(response2 ~ (1|item) + (1|participant) + Pair + Block + Orthography, data = datax) 
mod2 <- lmer(response2 ~ (1|item) + (1|participant) + Block + Pair*Orthography, data = datax)
emm <- emmeans(mod2, pairwise ~ Orthography | Pair, pbkrtest.limit = 5084)
emm
```

Plotting effects with emmeans:

```{r}
datax$Order <- datax$Pair
mod2 <- lmer(response2 ~ (1|item) + (1|participant) + Block + Order*Orthography, data = datax)
summary(mod2)

emm_options(pbkrtest.limit = 5084)
emmip(mod2, Orthography ~ Order, CIs = TRUE, col = "black",
      linearg = list(linetype = "solid", lwd = 1), dotarg = list(size = 4.5), CIarg = list(lwd = 0.5, alpha = 1)) + 
  ggplot2::theme_bw(base_size = 20) +
    ggplot2::xlab("Order")+
    ggplot2::ylab("Estimated marginal mean")
```

Now running cumulative link mixed models, given that *response* is really an ordinal variable.

```{r}
datax$response2 <- as.ordered(datax$response2)

mod1 <- clmm(response2 ~ (1|item) + (1|participant) + Block + Pair + Orthography, link = "logit", data = datax) 
mod2 <- clmm(response2 ~ (1|item) + (1|participant) + Block + Pair*Orthography, link = "logit", data = datax) 
anova(mod1, mod2)
summary(mod2)

#Getting the log odds and their confidence intervals.
exp(coef(mod2))
exp(confint(mod2))

#Plotting estimated marginal means on the latent scale with emmeans; this (sort of) creates a continuous response scale from the ordinal data so the plot looks identical to an emmeans plot for the equivalent lme model.

emmeans(mod2, pairwise ~ Orthography | Pair, mode= "latent")

emmip(mod2, Orthography ~ Pair, mode = "latent", CIs = TRUE, col = "black",
      linearg = list(linetype = "solid", lwd = 1), dotarg = list(size = 4.5), CIarg = list(lwd = 0.5, alpha = 1)) + 
  ggplot2::theme_bw(base_size = 20) +
    ggplot2::xlab("Pair")+
    ggplot2::ylab("Estimated marginal mean")

# Plotting predicted probabilities following method explained at https://marissabarlaz.github.io/portfolio/ols/ 

ggpredictions_mod2 = data.frame(ggpredict(mod2, terms = c("Orthography", "Pair"), type = "fe"))
ggpredictions_mod2 
ggpredictions_mod2$x = factor(ggpredictions_mod2$x)
levels(ggpredictions_mod2$x) = c("[+schwa]", "[-schwa]")
colnames(ggpredictions_mod2)[c(1, 6,7)] =c("Orthography", "Response", "Pair")
ggpredictions_mod2

# For some reason, if the colour of the CI lines is not set to Orthography, the dodge positioning of the points and lines becomes independent, so the points and whiskers don't align properly. To fix this, the colour must be set to Orthography, but the plot can still be made black/white by setting scale_colour_hue to l(uminance)=00.

ggpredictions_mod2$Orthography <- factor(ggpredictions_mod2$Orthography, levels = c("[-schwa]", "[+schwa]"))

ggplot(ggpredictions_mod2, aes(x = Response, y = predicted)) +
  geom_point(aes(shape = Orthography, color = Orthography), size = 3.5, position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = Orthography), position = position_dodge(width = 0.5), width = 0.01, alpha = 1) + 
  facet_wrap(~Pair) +
  scale_x_discrete(labels = c('slow','same','fast')) +
  xlab("")+
  ylab("Estimated probability") + 
  theme_bw(base_size = 18) +
  scale_colour_hue(l=00)

ggplot(ggpredictions_mod2, aes(x = Response, y = predicted, fill = Orthography)) + geom_bar(position = "dodge", stat = "identity") + facet_wrap(~Pair) + theme_bw() 

ggplot(ggpredictions_mod2, aes(x = Orthography, y = predicted, fill = Response)) + geom_bar(position = "fill", stat = "identity") + facet_wrap(~Pair) + theme_bw() 

```










