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


## Import data -----------------------------------------------------------------

data_morph=as.data.frame(read_excel("C:/Data/Thèse/Maroc_19_09/Manip_morphing/Raw_data.xlsx", sheet="Summary"))
data_morph$Stim_fact<-as.factor(data_morph$Stim)
data_morph$Croco<-as.factor(data_morph$Croco)
data_morph$HP<-as.factor(data_morph$HP)
data_morph$Ordre_stim_fact<-as.factor(data_morph$Ordre_stim)


# Score -----------------------------------------------------------------

data_morph$Score <- ordered(data_morph$Score, levels = c("0", "1", "2", "3", "4", "5"))

Model_median_morocco_linear <- brm(Score ~ Stim_fact + (1|Croco), data = data_morph, family = gaussian(),
                            iter = 5000, warmup = 500, chains = 3, cores = 3, seed = 1639)

summary(Model_median_morocco_linear)
plot(Model_median_morocco_linear, variable = c("b_Intercept[1]"))
pp_check(Model_median_morocco_linear) + theme_bw()
conditional_effects(Model_median_morocco_linear) 

# Plot prediction
newdata_morphing_morocco = data.frame(Stim_fact=levels(data_morph$Stim_fact))
fit_morphing_morocco = fitted(Model_median_morocco_linear, newdata = newdata_morphing_morocco, re_formula=NA, summary=T, robust=T)
pred_morphing_morocco = cbind(newdata_morphing_morocco, fit_morphing_morocco)
colnames(pred_morphing_morocco) = c('Stim', 'fit', 'se', 'lwr', 'upr')
head(pred_morphing_morocco)


# Full MCMC results
fit_morphing_full_morocco = fitted(Model_median_morocco_linear, newdata = newdata_morphing_morocco, re_formula=NA, summary=F)
colnames(fit_morphing_full_morocco) = newdata_morphing_morocco$Stim_fact
library(reshape2)
pred_morphing_full_morocco <- melt(fit_morphing_full_morocco, id.vars = 1, measure.vars = c("0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "1"))
colnames(pred_morphing_full_morocco) = c('Line', 'Stim', 'fit')

pred_morphing_morocco$Stim <- factor(as.factor(pred_morphing_morocco$Stim), levels = c("0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "1"))
pred_morphing_full_morocco$Stim <- factor(as.factor(pred_morphing_full_morocco$Stim), levels = c("0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "1"))

Morocco_morphing <- ggplot() +
  geom_violin(pred_morphing_full_morocco, mapping = aes(x=Stim, y=fit, fill=Stim), color="white", alpha=0.4, trim=TRUE)+
  geom_point(pred_morphing_morocco, mapping = aes(x=Stim, y=fit, color=Stim)) +
  geom_errorbar(pred_morphing_morocco, mapping = aes(x = Stim, ymin = lwr, ymax = upr, color=Stim), width = 0.3)+
  scale_color_manual(values=c("#009988", "#059C93", "#0BA09E", "#11A4AA", "#16A8B5", "#1CABC0", "#22AFCC", "#27B3D7", "#2DB7E2", "#33BBEE")) +
  scale_fill_manual(values=c("#009988", "#059C93", "#0BA09E", "#11A4AA", "#16A8B5", "#1CABC0", "#22AFCC", "#27B3D7", "#2DB7E2", "#33BBEE")) +
  theme_classic()+
  guides(fill = F, color=F) +
  labs(x="Morphing proportion", y="Behavioral score") + scale_x_discrete(labels=c("10%","20%","30%", "40%", "50%", "60%", "70%", "80%", "90%", "100%"))
Morocco_morphing

newdata_morocco_contrats = expand.grid(Stim_fact = levels(data_morph$Stim_fact))
contrast_morocco = as.data.frame(fitted (Model_median_morocco_linear, newdata_morocco_contrats, re_formula=NA, summary=F))
colnames(contrast_morocco) = newdata_morocco_contrats$Stim_fact
head(contrast_morocco)

boundary = contrast_morocco$`0.9` - contrast_morocco$`0.8`
soundgen:::reportCI(quantile(boundary, probs = c(.5, .025, .975)),2)
mean((contrast_morocco$`0.9` - contrast_morocco$`0.8`)>0)

croc_vs_frog_catego = ((contrast_morocco$`1` + contrast_morocco$`0.9`)/2) - ((contrast_morocco$`0.8` +
                                                                                contrast_morocco$`0.7` + contrast_morocco$`0.6` + contrast_morocco$`0.5` + contrast_morocco$`0.4`
                                                                              + contrast_morocco$`0.3` + contrast_morocco$`0.2` + contrast_morocco$`0.1`)/8)
quantile(croc_vs_frog_catego, probs = c(.5, .025, .975))
mean(croc_vs_frog_catego > 0)