---
title: "Analysis"
output: html_document
date: "2023-10-31"
editor_options: 
  chunk_output_type: console
---

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

##Assignment of stød in stress pattern groups

```{r}
library(tidyverse)
data = read.csv('Data_NonceWords.csv', sep=";")


#remove invalid datapoints:
data <- filter(data, Consistent_wprediction == "TRUE" | Consistent_wprediction == "FALSE" )

data <- mutate(data, speaker = factor(speaker))

#Add participant mean for participant results in the plot:
data2 <- data %>%
    group_by(speaker, Stress) %>%
  count(Consistent_wprediction, speaker) %>%
  mutate(prop_speaker = n / sum(n) *100) %>%
  filter(Consistent_wprediction == "TRUE")

#Add the mean manually when it was 0
data2 <- data2 %>% ungroup() %>% add_row(speaker = "1", Stress= "ante" , Consistent_wprediction =TRUE, n = 0 , prop_speaker = 0)
data2 <- data2 %>% ungroup() %>% add_row(speaker = "6", Stress= "ante" , Consistent_wprediction =TRUE, n = 0 , prop_speaker = 0)
data2 <- data2 %>% ungroup() %>% add_row(speaker = "5", Stress= "ante" , Consistent_wprediction =TRUE, n = 0 , prop_speaker = 0)


#plot the three stress pattern groups
data %>%
 count(Stress, Consistent_wprediction ) %>%
  group_by(Stress) %>%
  mutate(prop = n / sum(n)*100) %>%
  filter(Consistent_wprediction == "TRUE") %>%
  ggplot(aes(x=Stress, y= prop)) +
  geom_col(position=position_dodge()) + 
  geom_text(aes(label=round(prop)), position = position_dodge(.7), vjust = -.2, hjust =-.2)+
  geom_point(data=data2, aes(x=Stress, y=prop_speaker, fill=Consistent_wprediction), shape = 21)+
   scale_fill_manual( values = c("#619CFF")) +
  theme(legend.position="none")+
  labs(x="Stress Pattern Group", y="Percent Consistent with Prediction")+
  scale_x_discrete(labels=c("ante" = "Antepenultimate", "pen" = "Penultimate",
                              "ult" = "Ultimate"))

#ggsave("resultat1_true.png", width = 6, height = 6)
```

#Interspeaker differences

```{r}
#Get nonce words with stress on the ultimate syllable
data_ult <- filter(data, Stress == "ult")

#Get the mean for these nonce words for each speaker
speakers_ult <- data_ult %>%
  count(speaker, Consistent_wprediction) %>%
  group_by(speaker) %>%
  mutate(prop = round(n / sum(n) *100))

#To sort the participants from highest to lowest proportion responses that is "TRUE":
speakers_ult <- speakers_ult %>%
  mutate(prop_true = if_else(Consistent_wprediction == "TRUE", prop,0))

#Add column where the proportion responses that are "TRUE" is 0 to that column
speakers_ult <- speakers_ult %>% ungroup() %>% add_row(Consistent_wprediction =FALSE, speaker = "9", n = 0 , prop= 0, prop_true = 0)
speakers_ult <- speakers_ult %>% add_row(Consistent_wprediction =FALSE, speaker = "10", n = 0 , prop= 0, prop_true = 0)

#plot the proportion "TRUE" and proportion "FALSE" for each speaker as stables, show the precentages and number of responses for each:
ggplot(speakers_ult, aes(x = reorder(speaker, -prop_true), y=prop, fill = Consistent_wprediction)) +
  geom_col(position=position_dodge())+
  ylim(0,100)+
  geom_text(aes(label=round(prop)), position = position_dodge(.7), vjust = -.2, color ="black")+
  geom_text(aes(label="N="), position = position_dodge(.7), vjust = 2, hjust =1, size=3, color ="black")+
  geom_text(aes(label=n), position = position_dodge(.7), vjust = 2, size=3, hjust= -0.1, color ="black")+
  theme(legend.position="bottom")+
  labs(x="Speaker", y="Percent of responses")+
  labs(fill = "Consistent with prediction")+
  scale_fill_discrete(type = c("#f1a340", "#998ec3"), labels=c("No", "Yes"))+
  scale_x_discrete(labels=c("E101" = "1", "E102" = "2",
                              "E103" = "3", "E104"="4", "E105"="5", "E106"="6", "E107"="7", "E108"="8", "E109"="9", "E110"="10"))

#ggsave("Hoeg_Figure3.tiff", width = 12, height =4)
```

```{r}
#Get nonce words with stress on the antepenultimate syllable
data_ante <- filter(data, Stress == "ante")

#Get the mean for these nonce words for each speaker
speakers_ante <- data_ante %>%
  count(Consistent_wprediction, speaker) %>%
  group_by(speaker) %>%
  mutate(prop = round(n / sum(n) *100))

#To sort the participants from highest to lowest proportion responses that is "TRUE":
speakers_ante <- speakers_ante %>%
  mutate(prop_true = if_else(Consistent_wprediction == "TRUE", prop,0))

#Add column where the proportion responses that are "TRUE" is 0 to that column
speakers_ante <- speakers_ante %>% ungroup() %>% add_row(Consistent_wprediction =TRUE, speaker = "E110", n = 0 , prop= 0, prop_true = 0)
speakers_ante <- speakers_ante %>% add_row(Consistent_wprediction =TRUE, speaker = "1", n = 0 , prop= 0, prop_true = 0)
speakers_ante <- speakers_ante %>% add_row(Consistent_wprediction =TRUE, speaker = "6", n = 0 , prop= 0, prop_true = 0)
speakers_ante <- speakers_ante %>% add_row(Consistent_wprediction =TRUE, speaker = "5", n = 0 , prop= 0, prop_true = 0)

         
#plot the proportion "TRUE" and proportion "FALSE" for each speaker as stables, show the precentages and number of responses for each:
ggplot(speakers_ante, aes(x = reorder(speaker, -prop_true), y=prop, fill = Consistent_wprediction)) +
  geom_col(position=position_dodge())+
  ylim(0,100)+
  geom_text(aes(label=round(prop)), position = position_dodge(.7), vjust = -.2, color ="black")+
  geom_text(aes(label="N="), position = position_dodge(.7), vjust = 1.5, hjust =1, size=3, color ="black")+
  geom_text(aes(label=n), position = position_dodge(.7), vjust = 1.5, size=3, hjust= -0.1, color ="black")+
  theme(legend.position="bottom")+
  labs(x="Speaker", y="Percent of responses")+
  labs(fill = "Consistent with prediction")+
  scale_fill_discrete(type = c("#f1a340", "#998ec3"), labels=c("No", "Yes"))+
  scale_x_discrete(labels=c("E101" = "1", "E102" = "2",
                              "E103" = "3", "E104"="4", "E105"="5", "E106"="6", "E107"="7", "E108"="8", "E109"="9", "E110"="10"))

#ggsave("Hoeg_Figure4.tiff", width = 12, height =4)

```

```{r}
#Get nonce words with stress on the penultimate syllable
data_pen <- filter(data, Stress == "pen")

#Get the mean for these nonce words for each speaker
speakers_pen <- data_pen %>%
  count(speaker, Consistent_wprediction) %>%
  group_by(speaker) %>%
  mutate(prop = round(n / sum(n) *100))

#To sort the participants from highest to lowest proportion responses that is "TRUE":
speakers_pen <- speakers_pen %>%
  mutate(prop_true = if_else(Consistent_wprediction == "TRUE", prop,0))

#plot the proportion "TRUE" and proportion "FALSE" for each speaker as stables, show the precentages and number of responses for each:
ggplot(speakers_pen, aes(x = reorder(speaker, -prop_true), y=prop, fill = Consistent_wprediction)) +
  geom_col(position=position_dodge())+
  ylim(0,100)+
  geom_text(aes(label=round(prop)), position = position_dodge(.7), vjust = -.2, color ="black")+
  geom_text(aes(label="N="), position = position_dodge(.7), vjust = 2, hjust =1, size=3, color ="black")+
  geom_text(aes(label=n), position = position_dodge(.7), vjust = 2, size=3, hjust= -0.1, color ="black")+
  theme(legend.position="bottom")+
  labs(x="Speaker", y="Percent of responses")+
  labs(fill = "Consistent with prediction")+
  scale_fill_discrete(type = c("#f1a340", "#998ec3"), labels=c("No", "Yes"))+
  scale_x_discrete(labels=c("E101" = "1", "E102" = "2",
                              "E103" = "3", "E104"="4", "E105"="5", "E106"="6", "E107"="7", "E108"="8", "E109"="9", "E110"="10"))

#ggsave("Hoeg_Figure5.tiff", width = 12, height =4)
```


```{r}
#Remove speaker 1 to see how it affects the result:
data <- filter(data, speaker != "10")

#Med kun TRUE
data %>%
 count(Stress, Consistent_wprediction ) %>%
  group_by(Stress) %>%
  mutate(prop = n / sum(n)*100) %>%
  filter(Consistent_wprediction == "TRUE") %>%
  ggplot(aes(x=Stress, y= prop)) +
  geom_col(position=position_dodge()) + 
  geom_text(aes(label=round(prop)), position = position_dodge(.7), vjust = -.2, hjust =-.2)+
  geom_point(data=data2, aes(x=Stress, y=prop_speaker, fill=Consistent_wprediction), shape = 21)+
   scale_fill_manual( values = c("#619CFF")) +
  theme(legend.position="none")+
  labs(x="Stress Pattern Group", y="Percent Consistent with Prediction")+
  scale_x_discrete(labels=c("ante" = "Antepenultimate", "pen" = "Penultimate",
                              "ult" = "Ultimate"))

#ggsave("resultat1_true_speakerremoved.png", width = 6, height = 6)
```

#Item differences
```{r}
#Get nonce words with stress on the ultimate syllable after speaker 1 has been removed
data_ult <- filter(data, Stress == "ult")

#Get the number "TRUE" responses for each nonce word and the total number of responses.
words_ult <- data_ult %>%
  count(Nonceword, Consistent_wprediction) %>%
  group_by(Nonceword) 
words_ult <- words_ult %>%
  mutate(n_true = if_else(Consistent_wprediction == "TRUE", n,0))
words_ult_n <- data_ult %>%
  count(Nonceword) %>%
  mutate(n_total = n)
words_ult_n <- words_ult_n %>% select(Nonceword, n_total)
words_ult <- merge(words_ult, words_ult_n, by = "Nonceword")

#Order based on number of "TRUE" responses:
words_ult <- words_ult %>%
  arrange(-n_true)

#Separate the nonce words into two groups based on number of "TRUE" to split them into two plots:
groups = rep(1:2, each=56) #change each
words_ult$group = groups[1:nrow(words_ult)]

words_ult <- words_ult %>%
  arrange(Nonceword)

words_ult[c(3,7,11,21,23,27,30,32,36,44,68,70,75,77,79,80,88,93,102,104,106,108,117,119,125,127,129,134,135),6] <- 2

#Add "bigTU" manually as it had no valid datapoints:
words_ult[c(135),1] <- "bigTU"
words_ult[c(135),2] <- "FALSE"
words_ult[c(135),3] <- 0
words_ult[c(135),4] <- 0
words_ult[c(135),5] <- 0


words_ult <- words_ult %>%
  mutate(n_false = if_else(Consistent_wprediction == "FALSE", n,0))

#?
words_ult[c(135),7] <- 4

#Reorder based on number of "TRUE" and "FALSE":
words_ult <- words_ult %>%
  mutate(Nonceword=fct_reorder(Nonceword, -n_true, .fun='sum'))
words_ult <- words_ult %>%
  mutate(Nonceword=fct_reorder(Nonceword,n_false, .fun='sum'))

#Plot the first half of the nonce words:
ultplot1 <- ggplot(words_ult %>% filter(group =="1"), aes(x = Nonceword, y=n, fill = Consistent_wprediction)) +
  geom_col()+
  theme(axis.text.x=element_text(angle=98))+
  theme(text=element_text(size=18))+
  scale_fill_manual( values = c("#998ec3")) +
  labs(x=NULL, y="N of responses")+
  theme(legend.position="none")+
  coord_cartesian(clip = "off")+
  scale_y_continuous(breaks=seq(0,9,1))+
  geom_hline(yintercept = 1:9)

#Plot the second half of the nonce words:
ultplot2 <- ggplot(words_ult %>% filter(group =="2"), aes(x = Nonceword, y=n, fill = Consistent_wprediction)) +
  geom_col()+
  theme(axis.text.x=element_text(angle=100))+
  theme(text=element_text(size=18))+
  theme(legend.position="bottom")+
  labs(x="Nonce word", y="N of responses")+
  labs(fill = "Consistent with prediction")+
  scale_fill_discrete(type = c("#f1a340", "#998ec3"), labels=c("No", "Yes"))+
  coord_cartesian(clip = "off")+
  scale_y_continuous(breaks=seq(0,9,1))+
  geom_hline(yintercept = 1:9)

#Plot both together:
library(gridExtra)
library(cowplot)
grid.arrange(ultplot1, ultplot2, nrow=2)
plot_grid(ultplot1, ultplot2, align = "v", nrow = 2, rel_heights = c(3.2/7, 3.8/7))

#ggsave(file="Hoeg_Figure7.tiff", width = 20, height =12)
```


```{r}
#Get nonce words with stress on the antepenultimate syllable after speaker 1 has been removed
data_ante <- filter(data, Stress == "ante")

#Get the number "TRUE" responses for each nonce word and the total number of responses.
words_ante <- data_ante %>%
  count(Nonceword, Consistent_wprediction) %>%
  group_by(Nonceword) 
words_ante <- words_ante %>%
  mutate(n_true = if_else(Consistent_wprediction == "TRUE", n,0))
words_ante_n <- data_ante %>%
  count(Nonceword) %>%
  mutate(n_total = n)
words_ante_n <- words_ante_n %>% select(Nonceword, n_total)
words_ante <- merge(words_ante, words_ante_n, by = "Nonceword")

#Add "JUDtofe" manually as it had no valid datapoints:
words_ante <- words_ante %>% ungroup() %>% add_row(Nonceword = "JUDtofe", Consistent_wprediction= TRUE , n =0, n_true = 0, n_total = 0)


words_ante <- words_ante %>%
  mutate(n_false = if_else(Consistent_wprediction == "FALSE", n,0))

#to get "JUDtofe" in the rightmost position:
words_ante[107,6] <- 5.1


words_ante <- words_ante %>%
  mutate(Nonceword=fct_reorder(Nonceword, -n_total, .fun='sum'))

words_ante <- words_ante %>%
  mutate(Nonceword=fct_reorder(Nonceword, n_false, .fun='sum'))


#Plot the number of "TRUE" and "FALSE" for each nonce word, ordered from lowest to highest number of responses not consistent with the prediction:
ggplot(words_ante, aes(x =reorder(Nonceword,n_false, FUN=sum), y=n, fill = Consistent_wprediction))+
  geom_col()+
  theme(axis.text.x=element_text(angle=100))+
  theme(text=element_text(size=12))+
  theme(text=element_text(size=20))+
  theme(legend.position="bottom")+
  labs(x="Nonce word", y="N of responses")+
  labs(fill = "Consistent with prediction")+
  scale_fill_discrete(type = c("#f1a340", "#998ec3"), labels=c("No", "Yes"))+
  coord_cartesian(clip = "off")+
  scale_y_continuous(breaks=seq(0,9,1))+
  geom_hline(yintercept=1:5)
  

#ggsave(file="Hoeg_Figure8.tiff", width = 15, height =5)
```

```{r}
#Same precedure as for the nonce words with stress on the ultimate syllable.

data_pen <- filter(data, Stress == "pen")

words_pen <- data_pen %>%
  count(Consistent_wprediction, Nonceword) %>%
  group_by(Nonceword) 


words_pen <- words_pen %>%
  mutate(n_true = if_else(Consistent_wprediction == "TRUE", n,0))

words_pen_n <- data_pen %>%
  count(Nonceword) %>%
  mutate(n_total = n)

words_pen_n <- words_pen_n %>% select(Nonceword, n_total)

words_pen <- merge(words_pen, words_pen_n, by = "Nonceword")

words_pen <- words_pen %>%
  mutate(n_false = if_else(Consistent_wprediction == "FALSE", n,0))

words_pen <- words_pen %>%
  arrange(-n_true)

words_pen[c(1:29),7] <- 1
words_pen[c(31,33,35,36,37,38,39,40,41,42,43,44,46,48,49,50,51,52,53,54,56,58,60,61,67,68,71,72,79,80,85,89,90,91,96,102,110,118,123,127,128,129,131,132,133,137,138,139,142,143,145,146,147,148,149,150),7] <- 1

words_pen <- words_pen %>%
  arrange(V7)

words_pen[c(86:220),7] <- 2

#groups = rep(1:2, each=77)
#words_pen$group = groups[1:nrow(words_pen)]

words_pen <- words_pen %>%
  arrange(Nonceword)


words_pen <- words_pen %>% ungroup() %>% add_row(Nonceword = "SLIJril", Consistent_wprediction= TRUE , n =0, n_true = 0, n_total = 0, n_false = 0, V7 = 2)

words_pen <- words_pen %>% add_row(Nonceword = "GLÅMsuj", Consistent_wprediction= TRUE , n =0, n_true = 0, n_total = 0, n_false = 0, V7 = 2)

words_pen <- words_pen %>% add_row(Nonceword = "GURpev", Consistent_wprediction= TRUE , n =0, n_true = 0, n_total = 0, n_false = 0, V7 = 2)

words_pen <- words_pen %>% add_row(Nonceword = "HUboj", Consistent_wprediction= TRUE , n =0, n_true = 0, n_total = 0, n_false = 0, V7 = 2)

words_pen <- words_pen %>% add_row(Nonceword = "PRUNbir", Consistent_wprediction= TRUE , n =0, n_true = 0, n_total = 0, n_false = 0, V7 = 2)

words_pen <- words_pen %>% add_row(Nonceword = "STRÅNbir", Consistent_wprediction= TRUE , n =0, n_true = 0, n_total = 0, n_false = 0, V7 = 2)

words_pen <- words_pen %>% add_row(Nonceword = "TÅDraj", Consistent_wprediction= TRUE ,n =0, n_true = 0, n_total = 0, n_false = 0, V7 = 2)

words_pen <- words_pen %>%
  mutate(prop =  if_else(Consistent_wprediction == "TRUE", n_true / n_total * 100,0))

#words_pen[107,4] <- 0.0000001



#för att få ordningen att bli rätt i plot2. Den reorder som kommer sidst er den som "styrer".


words_pen <- words_pen %>%
  mutate(Nonceword=fct_reorder(Nonceword, -n_true, .fun='sum'))


penplot1 <- ggplot(words_pen %>% filter(V7 == "1"), aes(x = reorder(Nonceword,-n_true, FUN=sum), y=n, fill = Consistent_wprediction, )) +
  geom_col()+
  theme(axis.text.x=element_text(angle=98))+
  theme(text=element_text(size=20))+
  scale_fill_manual( values = c("#998ec3")) +
  labs(x=NULL, y="N of responses")+
  theme(legend.position="none")+
  coord_cartesian(clip = "off")+
  scale_y_continuous(breaks=seq(0,9,1))+
  geom_hline(yintercept = 1:9)

penplot2 <- ggplot(words_pen %>% filter(V7 == "2"), aes(x = reorder(Nonceword, -prop, FUN=sum), y=n, fill = Consistent_wprediction)) +
  geom_col()+
  theme(axis.text.x=element_text(angle=98))+
  theme(text=element_text(size=20))+
  theme(legend.position="bottom")+
  labs(x="Nonce word", y="N of responses")+
  labs(fill = "Consistent with prediction")+
  scale_fill_discrete(type = c("#f1a340", "#998ec3"), labels=c("No", "Yes"))+
  coord_cartesian(clip = "off")+
  scale_y_continuous(breaks=seq(0,9,1))+
  geom_hline(yintercept = 1:9)

library(gridExtra)
library(cowplot)
grid.arrange(penplot1, penplot2, nrow=2)
plot_grid(penplot1, penplot2, align = "v", nrow = 2, rel_heights = c(3.2/7, 3.8/7))

#ggsave(file="Hoeg_Figure9.tiff", width = 25, height =12)
```

#Danishness ratings
```{r}
#Prepare the data for the Danishness ratings:
data_Danishness = read.csv("Danishness_rating_data.csv", sep=";")

data_Danishness <- data_Danishness %>% select(Nonceword, Rating, Participant)

#Remove "always low rating" participants. Five other participants have been removed from the file earlier for the same reason.
data_Danishness <-  filter(data_Danishness, Participant != 22)
data_Danishness <-  filter(data_Danishness, Participant != 28)

#Remove all rows with NA:
data_Danishness <- data_Danishness %>% drop_na()

#Get the mean rating for each nonce word:
data_Danishness <- data_Danishness %>%
  group_by(Nonceword) %>%
  mutate(mean_rating = mean(Rating))

data_Danishness <- data_Danishness %>% select(Nonceword, mean_rating)

data_Danishness <- data_Danishness %>%
  group_by(Nonceword) %>%
  slice(1L)
```

Model the effect of Danishness rating of a nonce word have on response (with or without stød)
```{r}
#For nonce words with stress on the ultimate syllable
data_ult <- filter(data, Stress == "ult")

data_ult <- merge(data_ult,data_Danishness, by="Nonceword")

ggplot(data_ult, aes(Consistent_wprediction, mean_rating)) + geom_boxplot()


#converts "TRUE" to 1 and "FALSE" to 0:
data_ult$Consistent_wprediction <- as.integer(as.logical(data_ult$Consistent_wprediction))

danishness_ult_model <- glm(Consistent_wprediction ~ mean_rating, family=binomial(link = "logit"), data = data_ult)
summary(danishness_ult_model)


#predicted probability of a prediction-consistent response for lowest possible rating: 
#Based on the output obtained from the model when writing the paper. For some reason, the output for the oxytones is different now, but not the output for the paroxytones or the proparoxytones.
plogis(1.8363 + 0.5881 *1)
#0.9186691 (91% of the time)

plogis(1.8363 + 0.5881 *5)
#predicted probability of a prediction-consistent response for highest possible rating: 
#0.9916474 (99% of the time)
#thus, higher probability for a prediction-consistent response with a high Danishness rating
```


```{r}
#For nonce words with stress on the antepenultimate syllable:
data_ante <- filter(data, Stress == "ante")

data_ante <- merge(data_ante,data_Danishness, by="Nonceword")

ggplot(data_ante, aes(Consistent_wprediction, mean_rating)) + geom_boxplot()


#converts "TRUE" to 1 and "FALSE" to 0:
data_ante$Consistent_wprediction <- as.integer(as.logical(data_ante$Consistent_wprediction))

danishness_ante_model <- glm(Consistent_wprediction ~ mean_rating, family=binomial(link = "logit"), data = data_ante)
summary(danishness_ante_model)

#predicted probability of a prediction-consistent response for lowest possible rating: 
#estimate+estimate*lowest possible rating
plogis(-1.6428+0.4173*1)
#0.22697 (22%)

#predicted probability of a prediction-consistent response for highest possible rating: 
plogis(-1.6428+0.4173*5)
#0.6091403 (60%)
#thus, higher probability for a prediction-consistent response with a high Danishness rating, but not significant effect. 
#for et nonsensord med rating på 4, endast 50%
```


```{r}
#For nonce words with stress on the penultimate syllable:
data_pen <- filter(data, Stress == "pen")

data_pen <- merge(data_pen,data_Danishness, by="Nonceword")

ggplot(data_pen, aes(Consistent_wprediction, mean_rating)) + geom_boxplot()


#converts "TRUE" to 1 and "FALSE" to 0:
data_pen$Consistent_wprediction <- as.integer(as.logical(data_pen$Consistent_wprediction))

danishness_pen_model <- glm(Consistent_wprediction ~ mean_rating, family=binomial(link = "logit"), data = data_pen)
summary(danishness_pen_model)

#predicted probability of a prediction-consistent response for lowest possible rating: 
plogis(1.6568+0.1962*1)
#0.864479 (86%)

#predicted probability of a prediction-consistent response for highest possible rating: 
#0.9332551 (93%)

```