library(readxl)
library(ggplot2)
library(ggthemes)
library(brms)

labeling_complement = as.data.frame(read_excel("C:/Data/Thèse/Categorical_perception/Manip_labelling/Data_supplementary_morphing.xlsx", sheet="Raw"))
labeling_complement$Signal <- as.factor(labeling_complement$Signal)
labeling_complement$Croco <- as.factor(labeling_complement$Croco)


## Complementary experiments --------------------------------------------------- 

labeling_complement_test <- subset(labeling_complement, Signal!="CROC" & Signal!="FROG")
labeling_complement_test$Pourcentage <- factor(as.factor(labeling_complement_test$Pourcentage), levels = c("15", "20", "25", "30", "35", "40"))
labeling_complement_test$Type <- factor(as.factor(labeling_complement_test$Type), levels = c("full", "env"))

get_prior(Succes ~ Type + Pourcentage + (Pourcentage|Croco), data = labeling_complement_test, family = bernoulli("logit"))

library("psych")
Intercept =  rnorm(n = 1000, mean = 0, sd = 2)
hist(logistic(Intercept))
Effect_stim = rnorm(n= 1000, mean = 0, sd = 3)
hist(logistic (Intercept + Effect_stim) )

prior_comp = c(set_prior('normal(-2, 2)', class = 'Intercept'),
             set_prior('normal(-2, 2)', class = 'b', coef = "Pourcentage20"),
             set_prior('normal(0, 3)', class = 'b', coef = "Pourcentage25"),
             set_prior('normal(0, 3)', class = 'b', coef = "Pourcentage30"),
             set_prior('normal(0, 3)', class = 'b', coef = "Pourcentage35"),
             set_prior('normal(2, 2)', class = 'b', coef = "Pourcentage40"),
             set_prior('normal(0, 2)', class = 'b', coef = "Typeenv"))

labeling_complement_test$Pourcentage_num <- as.numeric(labeling_complement_test$Pourcentage)
Rep_morphing_comp <- brm(Succes ~ Type + Pourcentage + (Pourcentage|Croco), data = labeling_complement_test, family = bernoulli("logit"),
                         prior = prior_comp, iter=5000, warmup = 500, chains = 4, cores = 4, control = list(adapt_delta = 0.95), seed=1637) 

labeling_complement_test$Pourcentage_num <- as.numeric(as.character(labeling_complement_test$Pourcentage))
Rep_morphing_comp_sigmoid <- brm(Succes ~ Type + Pourcentage_num + (Pourcentage|Croco), data = labeling_complement_test, family = bernoulli("logit"),
                         iter=5000, warmup = 500, chains = 4, cores = 4, control = list(adapt_delta = 0.95), seed=1637) 

summary(Rep_morphing_comp)
pp_check(Rep_morphing_comp) + theme_bw()
conditional_effects(Rep_morphing_comp)


# Summary model predictions

newdata_morphing_comp =  expand.grid(Type=levels(labeling_complement_test$Type), Pourcentage=levels(labeling_complement_test$Pourcentage))
fit_morphing_comp = fitted(Rep_morphing_comp, newdata = newdata_morphing_comp, re_formula=NA, summary=T, robust=T)
colnames(fit_morphing_comp) = c('fit', 'se', 'lwr', 'upr')
Signal <- c("FULL_15", "ENV_15", "FULL_20", "ENV_20", "FULL_25", "ENV_25", "FULL_30", "ENV_30", "FULL_35", "ENV_35", "FULL_40", "ENV_40")
pred_morphing_comp = cbind(Signal, newdata_morphing_comp, fit_morphing_comp, Percent = c(15,15,20,20,25,25,30,30,35,35,40,40))

newdata_morphing_comp_sigmoid =  expand.grid(Type=levels(labeling_complement_test$Type), Pourcentage_num=seq(min(labeling_complement_test$Pourcentage_num), max(labeling_complement_test$Pourcentage_num), length.out = 100))
fit_morphing_comp_sigmoid = as.data.frame(fitted(Rep_morphing_comp_sigmoid, newdata = newdata_morphing_comp_sigmoid, re_formula=NA, summary=T, robust=T))
colnames(fit_morphing_comp_sigmoid) = c('fit', 'se', 'lwr', 'upr')
pred_morphing_comp_sigmoid = cbind(newdata_morphing_comp_sigmoid, fit_morphing_comp_sigmoid) 


# Full MCMC results
fit_morphing_comp_full = fitted(Rep_morphing_comp, newdata = newdata_morphing_comp, re_formula=NA, summary=F)
colnames(fit_morphing_comp_full) = c('FULL_15', 'ENV_15', 'FULL_20', 'ENV_20', 'FULL_25', 'ENV_25', 'FULL_30', 'ENV_30', 'FULL_35', 'ENV_35', 'FULL_40', 'ENV_40')
library(reshape2)
pred_morphing_comp_full <- melt(fit_morphing_comp_full, id.vars = 1, measure.vars = c('FULL_15', 'ENV_15', 'FULL_20', 'ENV_20', 'FULL_25', 'ENV_25', 'FULL_30', 'ENV_30', 'FULL_35', 'ENV_35', 'FULL_40', 'ENV_40'))
colnames(pred_morphing_comp_full) = c('Line', 'Signal', 'Fit')
pred_morphing_comp_full$Pourcentage <- pred_morphing_comp_full$Signal 
pred_morphing_comp_full$Pourcentage <- as.character(pred_morphing_comp_full$Pourcentage)
pred_morphing_comp_full$Pourcentage[pred_morphing_comp_full$Signal == "FULL_15" | pred_morphing_comp_full$Signal == "ENV_15"] <- "15"
pred_morphing_comp_full$Pourcentage[pred_morphing_comp_full$Signal == "FULL_20" | pred_morphing_comp_full$Signal == "ENV_20"] <- "20" 
pred_morphing_comp_full$Pourcentage[pred_morphing_comp_full$Signal == "FULL_25" | pred_morphing_comp_full$Signal == "ENV_25"] <- "25" 
pred_morphing_comp_full$Pourcentage[pred_morphing_comp_full$Signal == "FULL_30" | pred_morphing_comp_full$Signal == "ENV_30"] <- "30" 
pred_morphing_comp_full$Pourcentage[pred_morphing_comp_full$Signal == "FULL_35" | pred_morphing_comp_full$Signal == "ENV_35"] <- "35" 
pred_morphing_comp_full$Pourcentage[pred_morphing_comp_full$Signal == "FULL_40" | pred_morphing_comp_full$Signal == "ENV_40"] <- "40" 
pred_morphing_comp_full$Type <- pred_morphing_comp_full$Signal
pred_morphing_comp_full$Type <- as.character(pred_morphing_comp_full$Type)
pred_morphing_comp_full$Type[pred_morphing_comp_full$Signal == "FULL_15" | pred_morphing_comp_full$Signal == "FULL_20" | pred_morphing_comp_full$Signal == "FULL_25" | pred_morphing_comp_full$Signal == "FULL_30" | pred_morphing_comp_full$Signal == "FULL_35" | pred_morphing_comp_full$Signal == "FULL_40"] <- "full"
pred_morphing_comp_full$Type[pred_morphing_comp_full$Signal == "ENV_15" | pred_morphing_comp_full$Signal == "ENV_20" | pred_morphing_comp_full$Signal == "ENV_25" | pred_morphing_comp_full$Signal == "ENV_30" | pred_morphing_comp_full$Signal == "ENV_35" | pred_morphing_comp_full$Signal == "ENV_40"] <- "env" 

pred_morphing_comp_full$Type <- factor(as.factor(pred_morphing_comp_full$Type), levels = c("full", "env")) 
pred_morphing_comp$Type <- factor(as.factor(pred_morphing_comp$Type), levels = c("full", "env")) 


Plot_morphing_compl <- ggplot() +
  geom_violin(pred_morphing_comp_full, mapping = aes(x=Pourcentage, y=Fit, fill=Type, color=Type), alpha=0.4, color="white", position = position_dodge(width = 0.5), trim=TRUE)+
  geom_point(pred_morphing_comp, mapping = aes(x=Pourcentage, y=fit, color=Type), position = position_dodge(0.5)) +
  geom_errorbar(pred_morphing_comp, mapping = aes(x = Pourcentage, ymin = lwr, ymax = upr, color=Type), width = 0.2, position = position_dodge(0.5))+
  geom_smooth(pred_morphing_comp, mapping = aes(x=as.numeric(Pourcentage), y=fit, color=Type), method = "glm", method.args = list(family = "binomial"), se=FALSE, size=0.7) +
  #geom_line(pred_morphing_comp_sigmoid, mapping = aes(x=Pourcentage_num, y=fit, color=Type), size=0.7) +
  scale_color_manual(values=c("#00A9FF", "#F8766D"))+
  scale_fill_manual(values=c("#00A9FF", "#F8766D"), labels=c("whole", "spectral envelope"))+
  theme_classic()+
  guides(color = F, fill=F) +
  labs(x="Morphing proportion (%)", y="Response probability")
Plot_morphing_compl


# Detection threshold 
pred_morphing_comp$Pourcentage_num <- as.numeric(as.character(pred_morphing_comp$Pourcentage))
pred_morphing_comp_subfull <- subset(pred_morphing_comp, Type=="full")
pred_morphing_comp_subenv <- subset(pred_morphing_comp, Type=="env")
thresh_full = (max(pred_morphing_comp_subfull$fit)-min(pred_morphing_comp_subfull$fit))/2
thresh_env = (max(pred_morphing_comp_subenv$fit)-min(pred_morphing_comp_subenv$fit))/2

approx(pred_morphing_comp_subfull$fit, pred_morphing_comp_subfull$Pourcentage_num, xout = thresh_full)$y 
approx(pred_morphing_comp_subenv$fit, pred_morphing_comp_subenv$Pourcentage_num, xout = thresh_env)$y 


## Contrasts : credible interval (median of posterior distribution in %)

newdata = expand.grid(Type = levels(labeling_complement_test$Type), Pourcentage = levels(labeling_complement_test$Pourcentage))
contrast = as.data.frame(fitted (Rep_morphing_comp, newdata, re_formula=NA, summary=F))
colnames(contrast) = c("FULL_15", "ENV_15", "FULL_20", "ENV_20", "FULL_25", "ENV_25", "FULL_30", "ENV_30", "FULL_35", "ENV_35", "FULL_40", "ENV_40")

quantile(contrast$ENV_40 * 100, probs = c(.5, .025, .975))

f15_vs_e15 = contrast$ENV_40 - contrast$FULL_40
soundgen:::reportCI(quantile(f15_vs_e15 * 100, probs = c(.5, .025, .975)),1)
mean((contrast$ENV_40 - contrast$FULL_40)>0) * 100

ENV_30_vs_35 = contrast$ENV_35 - contrast$ENV_30
soundgen:::reportCI(quantile(ENV_30_vs_35 * 100, probs = c(.5, .025, .975)),1)
mean((contrast[, 6] - contrast[, 5])>0) * 100

croc_vs_frog_catego = ((contrast$FULL_40 + contrast$FULL_35)/2) - ((contrast$FULL_30 + contrast$FULL_25 + contrast$FULL_20 + contrast$FULL_15)/4)
soundgen:::reportCI(quantile(croc_vs_frog_catego * 100, probs = c(.5, .025, .975)),1)  
mean(croc_vs_frog_catego > 0) * 100

croc_vs_frog_catego = ((contrast$ENV_40 + contrast$ENV_35)/2) - ((contrast$ENV_30 + contrast$ENV_25 + contrast$ENV_20 + contrast$ENV_15)/4)
soundgen:::reportCI(quantile(croc_vs_frog_catego * 100, probs = c(.5, .025, .975)),1)  
mean(croc_vs_frog_catego > 0) * 100

croc_vs_frog_catego = ((contrast$FULL_40 + contrast$FULL_35 + contrast$FULL_30 + contrast$FULL_25 + contrast$FULL_20 + contrast$FULL_15)/6) - ((contrast$ENV_40 + contrast$ENV_35 + contrast$ENV_30 + contrast$ENV_25 + contrast$ENV_20 + contrast$ENV_15)/6)
soundgen:::reportCI(quantile(croc_vs_frog_catego * 100, probs = c(.5, .025, .975)),1)  
mean(croc_vs_frog_catego > 0) * 100




## Reaction time ---------------------------------------------------------------

compl_RT <- subset(labeling_complement_test, TpsReaction != "-1")

RT_compl <- brm(TpsReaction ~ Type + Pourcentage + (1|Croco), data = compl_RT, family = shifted_lognormal(),
                iter=5000, warmup=500, chains = 4, cores = 4, seed=1616, control = list(adapt_delta = 0.995)) #, max_treedepth = 15)) 

summary(RT_compl)
pp_check(RT_compl) + theme_bw() + xlim(0,100)
conditional_effects(RT_compl)

# Full MCMC results
fit_morphing_comp_full = as.data.frame(fitted(RT_compl, newdata = newdata_morphing_comp, re_formula=NA, summary=F))
colnames(fit_morphing_comp_full) = c('FULL_15', 'ENV_15', 'FULL_20', 'ENV_20', 'FULL_25', 'ENV_25', 'FULL_30', 'ENV_30', 'FULL_35', 'ENV_35', 'FULL_40', 'ENV_40')

soundgen:::reportCI(quantile(fit_morphing_comp_full$ENV_40, probs = c(.5, .025, .975)), 1)

croc_vs_frog_catego = ((fit_morphing_comp_full$ENV_40 + fit_morphing_comp_full$ENV_35)/2) - ((fit_morphing_comp_full$ENV_30 + fit_morphing_comp_full$ENV_25 + fit_morphing_comp_full$ENV_20 + fit_morphing_comp_full$ENV_15)/4)
soundgen:::reportCI(quantile(croc_vs_frog_catego, probs = c(.5, .025, .975)),1)  
mean(croc_vs_frog_catego > 0) * 100

croc_vs_frog_catego = ((fit_morphing_comp_full$FULL_40 + fit_morphing_comp_full$FULL_35)/2) - ((fit_morphing_comp_full$FULL_30 + fit_morphing_comp_full$FULL_25 + fit_morphing_comp_full$FULL_20 + fit_morphing_comp_full$FULL_15)/4)
soundgen:::reportCI(quantile(croc_vs_frog_catego, probs = c(.5, .025, .975)),1)  
mean(croc_vs_frog_catego > 0) * 100

croc_vs_frog_catego = ((fit_morphing_comp_full$FULL_40 + fit_morphing_comp_full$FULL_35 + fit_morphing_comp_full$FULL_30 + fit_morphing_comp_full$FULL_25 + fit_morphing_comp_full$FULL_20 + fit_morphing_comp_full$FULL_15)/6) - ((fit_morphing_comp_full$ENV_40 + fit_morphing_comp_full$ENV_35 + fit_morphing_comp_full$ENV_30 + fit_morphing_comp_full$ENV_25 + fit_morphing_comp_full$ENV_20 + fit_morphing_comp_full$ENV_15)/6)
soundgen:::reportCI(quantile(croc_vs_frog_catego, probs = c(.5, .025, .975)),1)  
mean(croc_vs_frog_catego > 0) * 100