---
title: "UsableAnalysis"
author: "Canaan Breiss"
date: "8/9/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
setwd("C:/Users/canaa/Google Drive/UCLA/Research and thoughts/VariationAndAdditivity/Analysis")
```

Load packages 

```{r}

#library(shinystan)

#library(sjstats,sjmisc)
library(lme4)
library(effects)
#library(strengejacke)
library(lmerTest)
library(gdata)
#library(ggplot2)
library(DescTools)
library(viridis)
library(car)
#library(dplyr)
#library(reshape2)
library(brms)
library(bayestestR)
library(tidyverse)
library(parameters)
```


Set up the dataframe

```{r}

######################################################
# ratings analysis
setwd("C:/Users/canaa/Google Drive/UCLA/Research and thoughts/VariationAndAdditivity/Analysis")
ratings <- read_csv("CombinedFinalDataset.csv")
colnames(ratings)[1] <- c("ObsNum")
ratings[7] <- NULL
ratings$Break <- as.factor(ratings$Break)
ratings$Conform <- as.factor(ratings$Conform)
ratings$TrialItem <- as.factor(ratings$TrialItem)
ratings$ExperimentName <- as.factor(ratings$ExperimentName)


ratings$Subject <- as.factor(ratings$Subject)
ratings$Rating <- as.numeric(ratings$Rating)

ratings <- drop.levels(subset(ratings, Rating <= 100))
ratings <- drop.levels(subset(ratings, Rating >= 0))

agg_ratings_means <- aggregate (Rating ~ Subject, ratings, mean)
agg_ratings_sd <- aggregate (Rating ~ Subject, ratings, sd)
names(agg_ratings_means)[2] <- "MeanRating"
names(agg_ratings_sd)[2] <- "SDofRating"
agg_rating <- cbind(agg_ratings_means, agg_ratings_sd[2])
big_df <- merge(ratings, agg_rating, by="Subject")
big_df <-cbind(big_df, ((big_df$Rating - big_df$MeanRating)/big_df$SDofRating))
names(big_df)[10] <- "Z"
big_df$Z <- as.numeric(big_df$Z)

ratings_addings <- read.csv("AddInsRatingFixed_NasalBackness.csv", sep = "\t",header = T)
ratings_addings[4] <- NULL
colnames(ratings_addings) <- c("Conform_real","Break_real","TrialItem","ViolatesBackness_real","ViolatesNasal_real","SameVowel_real","SameConsonants_real")

big_df <- merge(big_df, ratings_addings, by="TrialItem")
big_df$ViolatesBackness_real <- as.factor(big_df$ViolatesBackness_real)

big_df$ViolatesNasal_real <- as.factor(big_df$ViolatesNasal_real)
big_df$SameVowel_real <- as.factor(big_df$SameVowel_real)
big_df$SameConsonants_real <- as.factor(big_df$SameConsonants_real)
big_df$LogRating <- as.numeric(log10(big_df$Rating))
big_df$Subject <- as.factor((big_df$Subject))
big_df$TrialItem <- as.factor((big_df$TrialItem))
#big_df$Break_noneref <- relevel(big_df$Break_real, ref = "None" )
#big_df$Break_backnessref <- relevel(big_df$Break_real, ref = "Backness" )
#big_df$Break_nasalref <- relevel(big_df$Break_real, ref = "Nasal" )
#big_df$Break_backnesssibilantref <- relevel(big_df$Break_real, ref = "BacknessNasal" )

cleaned <- big_df 

cleaned$OrderedExperiments <- factor(cleaned$Experiment, levels=rev(levels(cleaned$Experiment)))
levels(cleaned$OrderedExperiments)


cleaned$NumExps  <- ifelse(cleaned$Experiment=="ExpQ100", "100",
                    ifelse(cleaned$Experiment=="ExpQ75", "75",
                    ifelse(cleaned$Experiment=="ExpQ50", "50", 
                    ifelse(cleaned$Experiment=="ExpQ63", "63",
                    ifelse(cleaned$Experiment=="ExpQ88", "88",NA)))))

cleaned$NumExps_regress  <- ifelse(cleaned$Experiment=="ExpQ100", "0",
                           ifelse(cleaned$Experiment=="ExpQ75", "25",
                                  ifelse(cleaned$Experiment=="ExpQ50", "50", 
                                         ifelse(cleaned$Experiment=="ExpQ63", "37",
                                                ifelse(cleaned$Experiment=="ExpQ88", "22",NA)))))

cleaned$NumExps <- as.numeric(cleaned$NumExps)
cleaned$NumExps_regress <- as.numeric(cleaned$NumExps_regress)

levels(cleaned$ViolatesNasal_real)

```
First the model for experiment 1

```{r}
m_ratings <- lmer(Rating~
                    #Break#*OrderedExperiments
                    #ViolatesNasal_real*Nasal_percent_correct+ViolatesBackness_real*Backness_percent_correct#*Experiment#+scale(TrialNumber)#+ViolatesBackness_real:Backness_percent_correct
                 +ViolatesNasal_real*ViolatesBackness_real*NumExps_regress
                  #+ ViolatesBackness_real*NumExps
                 +(1
                # +ViolatesNasal_real*ViolatesBackness_real
                |Subject)
                +(1|TrialItem)
                  ,data=cleaned);summary(m_ratings);plot(allEffects(m_ratings))
```

Now do the plot for experiment 1

```{r}
cleaned %>% 
  select(Subject,TrialItem,ExperimentName,Break_real,Rating) %>% 
  group_by(ExperimentName,Break_real) %>% 
  summarise(
    average = mean(Rating),
    sd = sd(Rating),
    se = se(Rating)
  ) %>% 
  mutate(Rating = average, 
         upper = average+se,
         lower = average-se,
         ExperimentsOrdered = factor(ExperimentName, levels=rev(levels(cleaned$ExperimentName)))) %>% 
ggplot(aes(y = Rating, ymin = lower, ymax=upper, x = ExperimentsOrdered, color = Break_real, group = Break_real))+
  geom_pointrange(size=0.25)+
  geom_line()+
  ylab("Rating\n (mean and standard error)")+
  xlab("Condition\n (Fraction exceptional forms for each phonotactic)")+
  scale_x_discrete(labels=c("A\n (0/32)","B\n (2/32)","C\n (4/32)","D\n (6/32)","E\n (8/32)"))+
  scale_color_viridis(option = "viridis", discrete = T, labels = c("Violates backness harmony","Violates both phonotactics","Violates nasal harmony","Violates neither phonotactic"), name = "Violation profile")+
  theme_bw()+
  ylim(30,75)+
  theme(legend.position = "right")


ggsave("Exp1.jpeg",height = 4, width=6)

```

Now try doing the faceted plots - full facets

```{r}

hospital_names <- list(
  'ExpQ100'="A\n (0/32)",
  'ExpQ88'="B\n (2/32)",
  'ExpQ75'="C\n (4/32)",
  'ExpQ63'="D\n (6/32)",
  "ExpQ50" = "E\n (8/32)"
)

hospital_labeller <- function(variable,value){
  return(hospital_names[value])
}
v <- cleaned$ExperimentName
names(v) <- c("A\n (0/32)","B\n (2/32)","C\n (4/32)","D\n (6/32)","E\n (8/32)")

l <- data.frame(c("A\n (0/32)","B\n (2/32)","C\n (4/32)","D\n (6/32)","E\n (8/32)"))

d <- cleaned %>% 
  select(Subject,TrialItem,ExperimentName,Rating,ViolatesBackness_real,ViolatesNasal_real) %>% 
  group_by(ExperimentName,ViolatesBackness_real,ViolatesNasal_real) %>% 
  summarise(
    average = mean(Rating),
    sd = sd(Rating),
    se = standard_error(Rating)
  ) %>% 
  mutate(Rating = average, 
         upper = average+se,
         lower = average-se,
         ExperimentsOrdered = factor(ExperimentName, levels=rev(levels(cleaned$ExperimentName))))# %>% 

v <- d$ExperimentsOrdered
names(v) <- c("A\n (0/32)","B\n (2/32)","C\n (4/32)","D\n (6/32)","E\n (8/32)")

d %>% 
ggplot(aes(y = Rating, ymin = lower, ymax=upper, x = ViolatesBackness_real, color = ViolatesNasal_real, group = ViolatesNasal_real))+
  geom_errorbar(size=0.25)+
  geom_line()+
  ylab("Rating\n (mean and standard error)")+
  xlab("Violation of backness harmony")+
  scale_x_discrete(labels=c("No", "Yes"))+
  scale_color_viridis(option = "viridis", discrete = T, labels = c("No","Yes"), name = "Violation of nasal harmony", begin = .3, end=0.7)+
  theme_bw()+
  ylim(30,75)+
  theme(legend.position = "bottom")+
  facet_grid(cols = vars(ExperimentsOrdered), 
             labeller = hospital_labeller)
        
ggsave("Exp1_Interactionplot_Full.jpeg", height = 5, width = 8)     
```





Nasal advantage

```{r}

v_50 <- read_csv("./VerifyTrials/usable_trials_lexdec_nasalback_distributional_verify_Q50.csv")
v_63 <- read_csv("./VerifyTrials/usable_trials_lexdec_nasalback_distributional_verify_Q63.csv")
v_75 <- read_csv("./VerifyTrials/usable_trials_lexdec_nasalback_distributional_verify_Q75.csv")
v_88 <- read_csv("./VerifyTrials/usable_trials_lexdec_nasalback_distributional_verify_Q88.csv")
v_100 <- read_csv("./VerifyTrials/usable_trials_lexdec_nasalback_distributional_verify_Q100.csv")

v_all <- rbind(v_100,v_88,v_75,v_63,v_50)

just_v_2_subjects <- v_all %>% 
  filter(Phase == "verify_2") %>% 
  dplyr::select(Subject) %>% 
  distinct()


just_v_2_trials <- v_all %>% 
  filter(Phase == "verify_2") 

just_v_1 <- v_all %>% 
  filter(!(Subject %in% list(just_v_2_subjects$Subject))) 
  

test <- v_all %>% 
  
  mutate(Correct_numeric = if_else(Correct == "right",1,0)) %>% 
  group_by(ExperimentName,Subject,Break) %>% 
  mutate(TotalCorrect = sum(Correct_numeric)) %>% 
  ungroup() %>% 
  dplyr::select(TotalCorrect,Subject,Break,ExperimentName) %>% 
  distinct() %>% 
  filter(TotalCorrect <9) #%>%
 


difference <- test %>% 
  pivot_wider(names_from = Break, values_from = TotalCorrect) %>%
  mutate(Subject = paste(Subject, ExperimentName, sep = "_")) %>% 
  mutate(Difference = Nasal-Backness) %>% 
  mutate(ExperimentName = as.factor(ExperimentName))# positive is more nasals correct

ggplot(difference, aes(y = Difference, x = ExperimentName, color = ExperimentName))+
  geom_jitter(alpha=.5, height=0.2, width=0.2)+
  ylab("Nasal advantage in verification phase")+
  ylim(-3.9,3.9)+
  xlab("Condition\n (Fraction exceptional forms for each phonotactic)")+
  scale_x_discrete(labels=c("A\n (0/32)","B\n (2/32)","C\n (4/32)","D\n (6/32)","E\n (8/32)"))+
  scale_color_viridis(option = "plasma", discrete = T)+
  theme_bw()+
  theme(legend.position = "none")

ggsave("NasalAdvantage.jpeg",height = 4, width=5)


```


```{r}
# check with regression
m_nasadvantage <- lm(Difference ~ PercentConforming,
                       data = difference %>% 
                       mutate(PercentConforming = str_replace(ExperimentName, "ExpQ",""),
                              PercentConforming = as.numeric(PercentConforming)))
summary(m_nasadvantage)


# This is to check to see if the effect still holds even if we restrict to people who were right in the middle - no underlearning likely.
difference_between_minus_one_and_one <- difference %>% 
  #filter(Difference == 0) %>% 
  filter((Difference > -2),(Difference < 2)) %>% 
  dplyr::select(Subject) %>% 
  distinct()
data_just_middles <- 
  subset(cleaned, Subject %in% difference_between_minus_one_and_one$Subject)#% 
  

m_just_middles <- lmer(Rating~
                   #Break#*OrderedExperiments
                   #ViolatesNasal_real*Nasal_percent_correct+ViolatesBackness_real*Backness_percent_correct#*Experiment#+scale(TrialNumber)#+ViolatesBackness_real:Backness_percent_correct
                   +ViolatesNasal_real*NumExps_regress*ViolatesBackness_real
                 
                 +(1
                   # +ViolatesNasal_real*ViolatesBackness_real
                   |Subject)
                 +(1|TrialItem)
                 ,data=data_just_middles);summary(m_just_middles);plot(allEffects(m_just_middles))

```


Now compare with exp 2

```{r}
ratings_A <- read_csv("usable_ratings_lexdec_nasalback_distributional_Q50A.csv")
ratings_B <- read_csv("usable_ratings_lexdec_nasalback_distributional_Q50B.csv")
ratings_A$ExperimentName <- "ExpQ50A"

ratings_combined <- rbind(ratings_A,ratings_B)
#colnames(ratings)[1] <- c("ObsNum")
#ratings[7] <- NULL
ratings_combined$Break <- as.factor(ratings_combined$Break)
ratings_combined$Conform <- as.factor(ratings_combined$Conform)
ratings_combined$TrialItem <- as.factor(ratings_combined$TrialItem)
ratings_combined$ExperimentName <- as.factor(ratings_combined$ExperimentName)


ratings_combined$Subject <- as.factor(ratings_combined$Subject)
ratings_combined$Rating <- as.numeric(ratings_combined$Rating)

ratings_combined <- drop.levels(subset(ratings_combined, Rating <= 100))
ratings_combined <- drop.levels(subset(ratings_combined, Rating >= 0))

agg_ratings_means_combined <- aggregate (Rating ~ Subject, ratings_combined, mean)
agg_ratings_sd_combined <- aggregate (Rating ~ Subject, ratings_combined, sd)
names(agg_ratings_means_combined)[2] <- "MeanRating"
names(agg_ratings_sd_combined)[2] <- "SDofRating"
agg_rating_combined <- cbind(agg_ratings_means_combined, agg_ratings_sd_combined[2])
big_df_combined <- merge(ratings_combined, agg_rating_combined, by="Subject")
big_df_combined <-cbind(big_df_combined, ((big_df_combined$Rating - big_df_combined$MeanRating)/big_df_combined$SDofRating))
names(big_df_combined)[10] <- "Z"
big_df_combined$Z <- as.numeric(big_df_combined$Z)

ratings_addings <- read.csv("AddInsRatingFixed_NasalBackness.csv", sep = "\t",header = T)
ratings_addings[4] <- NULL
colnames(ratings_addings) <- c("Conform_real","Break_real","TrialItem","ViolatesBackness_real","ViolatesNasal_real","SameVowel_real","SameConsonants_real")

big_df_combined <- merge(big_df_combined, ratings_addings, by="TrialItem")
big_df_combined$ViolatesBackness_real <- as.factor(big_df_combined$ViolatesBackness_real)

big_df_combined$ViolatesNasal_real <- as.factor(big_df_combined$ViolatesNasal_real)
big_df_combined$SameVowel_real <- as.factor(big_df_combined$SameVowel_real)
big_df_combined$SameConsonants_real <- as.factor(big_df_combined$SameConsonants_real)
big_df_combined$LogRating <- as.numeric(log10(big_df_combined$Rating))
big_df_combined$Subject <- as.factor((big_df_combined$Subject))
big_df_combined$TrialItem <- as.factor((big_df_combined$TrialItem))
#big_df_combined$Break_noneref <- relevel(big_df_combined$Break_real, ref = "None" )
#big_df_combined$Break_backnessref <- relevel(big_df_combined$Break_real, ref = "Backness" )
#big_df_combined$Break_nasalref <- relevel(big_df_combined$Break_real, ref = "Nasal" )
#big_df_combined$Break_backnesssibilantref <- relevel(big_df_combined$Break_real, ref = #"BacknessNasal" )

cleaned_combined <- big_df_combined

cleaned_combined$OrderedExperiments <- factor(cleaned_combined$Experiment, levels=rev(levels(cleaned_combined$Experiment)))
levels(cleaned_combined$OrderedExperiments)
```


Bayesian model


```{r}

b_ratings <- brm(Rating~
                    +ViolatesNasal_real*ViolatesBackness_real*OrderedExperiments
                  +(1
                     +ViolatesNasal_real*ViolatesBackness_real
                    |Subject)
                  +(1+OrderedExperiments|TrialItem),
                 data=cleaned_combined %>% 
                   select(
                     ViolatesBackness_real,ViolatesNasal_real,OrderedExperiments,Subject,TrialItem,Rating
                     ),
                 cores = 4, 
                 chains = 4, 
                 warmup = 1000, 
                 iter = 10000,
                 save_all_pars = T,
                 file = "bayesian_model_exp_2",
                 save_model = "bayesian_exp2_model_stan")
summary(b_ratings)
conditional_effects(b_ratings)
p_direction(b_ratings)

```
Now plot exp 2
```{r}

cleaned_combined %>% 
  select(Subject,TrialItem,ExperimentName,Break_real,Rating) %>% 
  group_by(ExperimentName,Break_real) %>% 
  summarise(
    average = mean(Rating),
    sd = sd(Rating),
    se = se(Rating)
  ) %>% 
  mutate(Rating = average, 
         upper = average+se,
         lower = average-se,
         ExperimentsOrdered = factor(ExperimentName, levels=(levels(cleaned_combined$ExperimentName)))) %>% 
ggplot(aes(y = Rating, ymin = lower, ymax=upper, x = ExperimentsOrdered, color = Break_real))+
  geom_pointrange(size=0.25)+
  ylab("Rating\n (mean and standard error)")+
  xlab("Experiment")+
  scale_x_discrete(labels=c("Experiment 1, Condition E","Experiment 2"))+
  scale_color_viridis(option = "viridis", discrete = T, labels = c("Violates backness harmony","Violates both phonotactics","Violates nasal harmony","Violates neither phonotactic"), name = "Violation profile")+
  ylim(30,75)+
  theme_bw()+
  theme(legend.position = "right")


ggsave("Exp2.jpeg",height = 4, width=6)

```

