---
title: "Noemie memoire"
author: "Madarasz Noemie"
date: "11/30/2022"
output: html_document
---

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

##Memoire de Noemie Madarasz sur la domestication
#Les poissons utilisés sont 42 poissons de la génération F0 et 28 poissons de la génération F5 de la lignée EPP de Kryptolebias marmoratus
#Ces poissons ont été soumis à un test comportemental mesurant l'audace avec deux répliquas séparés de 2j (shelter test)
#Ce test consiste en une arène ouverte accolée à un abri fermé dans lequel le poisson est positionné en début de test. On délimite dans l'arène une zone interne

#Le jeu de données utilisé (Memoire_EPP.csv) comprend : 
# - Fish ID 
# - Generation (F0 ou F5)
# - Replicate (1 ou 2)
# - Weight du poisson (g) après le deuxième répliqua
# - Length du poisson (cm) après le deuxième répliqua
# - Total time in open field (s) = TTA = Le temps total que le poisson a passé au sein de l'arène ouverte (zone interne + zone externe)
# - Latency to first entry in open field (s) = Latarena = La latence pour que le poisson sorte du shelter pour la première fois pour aller dans l'arène
# - Total time in internal zone (s) = TTI = Le temps total que le poisson a passé dans la zone interne délimitée dans l'arène
# - Total time in external zone (s) = TTE = Le temps total que le poisson a passé dans la zone externe, c'est à dire tout le reste de l'arène autre que la zone interne
# - Latency to first entry in internal zone (s) = Latint = La latence pour que le poisson entre dans la zone interne pour la première fois
# - Total distance moved (cm) = TDM = La distance totale parcourue dans le poisson dans l'arène
# - Total time in shelter (s)" = TTS = Le temps total que le poisson a passé dans l'abri (donc la différence du temps total du test - le temps total dans toute l'arène 

```{r}

library(sjPlot)
library(ggpubr)
library(psych)
library(Hmisc)
library(ade4)
library(vegan)
library(lme4)
library(optimx)
library(car)
library(AICcmodavg)
library(VGAMdata)
library(fitdistrplus)
library(ROCR)
library(MCMCglmm)
library(car)
library(readr)
library(RVAideMemoire)
library(tidyverse)
library(data.table)
library(lubridate)
library(lme4)
library(lmerTest)
library(rptR)
library(Hmisc)
library(bestNormalize)
library(knitr)
library(corrplot)
library(xtable)
library(PerformanceAnalytics)
library(kableExtra) 
library(factoextra) 
library(FactoMineR) 
library(ggplot2)
library(png) 
library(AlphaPart)
library(ggrepel)
library(cowplot)
library(ggspatial)
library(googleway)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(dplyr)
library("ade4")
library("factoextra")
library(tidyr)
library(pwr)
library(readxl)

```

============================================================================
                           Importation des données 
============================================================================

```{r}
ls()
rm(list = ls())

EPP <- read_excel("Memoire_EPP.xlsx", sheet="Memoire_EPP")
summary(EPP)

#renommer les variables
EPP <- rename(EPP, c("ID"="Fish ID","Weight"="Weight (g)","Length"="Length (cm)","TDM"="Total distance moved (cm)","TTA"="Total time in open field (s)","Latarena"="Latency to first entry in open field (s)","TTI"="Total time in internal zone (s)", "Latint"="Latency to first entry in internal zone (s)","TTS"="Total time in shelter (s)"))
View(EPP)
str(EPP)

#Transformer les variables non continues en facteur 
EPP$`ID`=as.factor(EPP$`ID`)
EPP$`Generation`=as.factor(EPP$`Generation`)
EPP$`Replicate`=as.factor(EPP$`Replicate`)

summary(EPP)
str(EPP)

## nouvelles variables BMI et RTDM
EPP <- mutate(EPP, RTDM=TDM/Length/TTA) #TDM est divisé par le TTA pour éviter que l'audace n'interfère avec l'activité car la distance parcourue est proportionnelle au tems passé dans l'arène 

EPP$IZ <- EPP$TTI/EPP$TTA #IZ = fraction du temps passé dans la zone interne (TTI/TTA)

EPP <- mutate(EPP, BMI=Weight/Length^2)

View(EPP)

EPP <- EPP[,c("Generation", "ID", "Replicate", "Weight", "Length", "BMI", "TTA", "Latarena", "Latint", "RTDM","IZ","TTI")] #subset only the columns of interest 

attach(EPP)

```

#Visualisation des variables de poids et de taille + BMI

```{r}
df <- EPP[1:70,1:6] 
Weight <- ddply(df, c("Generation"), summarize, Mean = mean(Weight), SD = sd(Weight))
Length <- ddply(df, c("Generation"), summarize, Mean = mean(Length), SD = sd(Length))
BMI <- ddply(EPP, c("Generation"), summarize, Mean = mean(BMI), SD = sd(BMI))

ggplot(BMI) +
  aes(x = Generation, y = Mean, fill=Generation) + 
  scale_fill_manual(breaks = c("1", "2"), values=c("darkgrey","black"))+
  geom_bar(position=position_dodge(0.65), colour="black", stat = "identity", width=0.55)  +
  theme_bw() +
  ylab("BMI") +
  theme(aspect.ratio = 1/1)+
  scale_y_continuous(limits=c(0, 0.06))+
  geom_errorbar(position=position_dodge(0.65), aes(ymin =Mean , ymax = Mean + SD), width = 0.55)


```

#============================================================================
#                  Transfo données et graphiques exploration des données
#============================================================================

```{r}

vars <- read_excel("Memoire_EPP.xlsx", sheet="variables")

EPP2 <- EPP %>% gather(variable, value, -Replicate, -Generation, -ID) %>%
  mutate(Vars = factor(variable, 
                       levels=vars$variable,
                       labels=vars$nom)) %>%
  merge( (vars %>% select(variable)))

EPP2


ggplot(EPP2, aes(y=value, colour=Generation))+
  facet_wrap(~Vars, scales="free_y")+
  geom_point(alpha=0.5, position=position_jitter(width=0.1), aes(x=Replicate))+
  theme_bw()+
  stat_summary(geom="line", fun="mean", aes(group=variable,x=Replicate))+
  stat_summary(geom="point", fun="mean", aes(group=variable,x=Replicate))


ggplot(EPP2, aes(y=value, colour=Replicate))+
  facet_wrap(~Vars, scales="free_y")+
  geom_point(alpha=0.5, position=position_jitter(width=0.1), aes(x=Generation))+
  theme_bw()+
  stat_summary(geom="line", fun="mean", aes(group=variable,x=Generation))+
  stat_summary(geom="point", fun="mean", aes(group=variable,x=Generation))
```

#============================================================================
#                                   CORRELATIONS
#============================================================================

```{r}
EPP
# Sélection des colonnes avec nos variables réponses et explicatives. Spearman = non gaussian, Pearson = gaussian
correlation <- rcorr(as.matrix(EPP[,4:11]), type = "spearman") 

# Résumé de la matrice : 3 listes : r = coef de correlation ; n = nombre de valeurs ; P = p valeur des corrélations
str(correlation) 

corrplot(correlation$r, type="upper", order="hclust", p.mat=correlation$p, sig.level=0.01, insig="blank", tl.cex=0.7 ,tl.col="black",tl.srt=45)

#+ c'est bleu + c'est corrélé positivement
#Length et Weight et BMI

heatmap(correlation$r, symm=T)

# Sélection du tableau avec les coef de corrélations (r) et du tableau avec les p valeurs (P)
r=data.frame(correlation$r) 
p=data.frame(correlation$P)

chart.Correlation(as.matrix(EPP[,4:11]), histogram = T,pch=19)

# flattenCorrMatrix function
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut])
}

flat = flattenCorrMatrix(correlation$r, correlation$P)
flat

# Exportation des tableaux
write_csv(r,'correlationmatrix.csv')
write_csv(p,'pvalcorrelationmatrix.csv')
write_csv(flat, 'matrix.csv')

```


#===============================================================================
#                           ACP sur les variables comp
#===============================================================================


## ---------------- GENERATIONS F0 ET F5 ; R1 et R2 -------------

```{r}
EPP <- na.omit(EPP) 
EPP
acp <- dudi.pca(EPP[,7:11], scale=TRUE,center=TRUE,nf = 3) #nf = nombre de composantes dans l'ACP
3
acp

#eigenvalues = valeurs propres
acp$eig
summary(acp) #Ax1 58.543  Ax2 19.281 de variabilité représentée par l'ACP 

barplot(acp$eig, xlab="Axe", ylab="Valeur propre") #Montre en barplot les valeurs propres
fviz_pca_biplot(acp)

screeplot <- fviz_eig(acp)
screeplot

ACP <- fviz_pca_biplot(acp, axes = c(1,2), geom.ind = "point", fill.ind = EPP$Generation,  col.ind = "black", pointshape = 21, pointsize = 2, palette = c("#999999", "#E69F00", "#56B4E9"), addEllipses = TRUE, col.var = "contrib", gradient.cols=c("red", "gold", "forestgreen"), select.var=list(contrib=10), repel = T, legend.title = list(fill = "Generation", color = "Contribution"))
ACP #très belle figure

#loadings
acp$c1
#correlations between variables and PCs (principal components)
acp$co
compo=acp$co
compo
str(compo)
write_csv(data.frame(compo),'compo.csv') #exporter

# REPRESENTATION 

fviz_pca_var(acp, axes=c(1,2), 
             title = "Cercle des correlations selon le plan 1-2",
             repel = T)

scatter(acp)#Représentation de tout (valeurs, composantes, et eigenvalues)

#to show contribution of each variable. Color "contrib" = Color by the contribution. repel = T to Avoid text overlapping
varcont <- fviz_pca_var(acp, col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),repel = T,)
varcont

#to show variables with the 3 highest contributions
varcont2 <- fviz_pca_var(acp,
                         col.var = "contrib", 
                         gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                         select.var=list(contrib=3),
                         repel = T,)
varcont2

#Contribution axe 1
contribution <- fviz_contrib(acp, choice = "var", axes = 1)
contribution

#Contribution axe 2
contribution2 <- fviz_contrib(acp, choice = "var", axes = 2)
contribution2

#Contribution axe 3 
contribution3 <- fviz_contrib(acp, choice = "var", axes = 3)
contribution3

ggdraw() + draw_plot(ACP, 0, .5, 1, .5) + draw_plot(contribution, 0, 0, .5, .5) + draw_plot(contribution2, .5, 0, .5, .5) + draw_plot_label(c("A", "B", "C"), c(0, 0, 0.5), c(1, 0.5, 0.5), size = 15)

#la dimension 1 est très forte
#sélection de TTA, IZ, RTDM et Latarena
#la variabilité des F0 est bcp plus grande que les F5




#Visualiser l'effet de Replicate VS Generation
#uniquement sur les variables comportementales (pas taille et poids)
acp2 <- dudi.pca(EPP[,7:11], scale=TRUE,center=TRUE,nf = 3) #nf = nombre de composantes dans l'ACP
3
acp2
s.class(acp2$li,EPP$Replicate,col=c("royalblue", "red3"),xax=2,yax=3,cel=0,sub="Replicate")
s.class(acp2$li,EPP$Generation,col=c("royalblue", "red3","green", "orange"),xax=2,yax=3,cel=0,sub="Generation")

#--> effet plus important de generation


```


#===============================================================================
#                           LINEAR MODELS
#===============================================================================
#On a choisi d'etudier trois parametres

#1) Si le poisson est audacieux et va volontairement et rapidement quitter le shelter pour explorer l'arene (Latency to first entry in arena)
#2) Si l'animal prend des risques et rest longtemps dans l'arene (Total time in arena)
#3) Si l'animal est anxieux et va avoir tendance à longer les parois et ne pas explorer le centre de l'arene, appelé phénomène de thigmotaxisme (IZ ; en proportion et independemment de l'audace)
#4) L'activité des poissons (RTDM ; independemment de l'audace)


# 1. ----------------- Latency to first entry in arena -----------------------------------------------

retirer weight et BMI (covarie avec length)

```{r}
#récupérer les lignes qui ont été supprimées à cause des NAs
EPP <- read_excel("Memoire_EPP.xlsx", sheet="Memoire_EPP")
EPP <- rename(EPP, c("ID"="Fish ID","Weight"="Weight (g)","Length"="Length (cm)","TDM"="Total distance moved (cm)","TTA"="Total time in open field (s)","Latarena"="Latency to first entry in open field (s)","TTI"="Total time in internal zone (s)", "Latint"="Latency to first entry in internal zone (s)","TTS"="Total time in shelter (s)"))
EPP$`ID`=as.factor(EPP$`ID`)
EPP$`Generation`=as.factor(EPP$`Generation`)
EPP$`Replicate`=as.factor(EPP$`Replicate`)
EPP <- mutate(EPP, RTDM=TDM/Length/TTA) 
EPP <- mutate(EPP, IZ=TTI/TTA)
EPP <- mutate(EPP, BMI=Weight/Length^2)

test1 <- EPP[,c("Generation", "ID", "Replicate", "Length", "Latarena")] #subset only the columns of interest
test1 <- na.omit(test1) #delete the rows with NA in all the columns
test1
summary(test1)
head(test1)

```


```{r}
test1
# search for extreme values
dotchart(test1$Latarena, main="Latency to first entry in arena", group=test1$Generation)

hist(test1$Latarena)
#visualiser 
#Pas normal, homogene pas suivie
#Pour remédier à ça : on transforme nos données

hist(test1$Length)
shapiro.test(test1$Length) #significatif, pas ok

bestNormalize(test1$Latarena) # prendre Yeo-J
bestNormalize(test1$Length) # > orderNorm > prendre Yeo-J

LA <- yeojohnson(test1$Latarena)
x1 <- predict(LA)
test1$Latarena <- x1
hist(test1$Latarena)
shapiro.test(test1$Latarena) # significatif mais moins -> ok, ça passe

TA <- yeojohnson(test1$Length)
x2 <- predict(TA)
test1$Length <- x2
hist(test1$Length)
shapiro.test(test1$Length) # significatif mais moins -> ok; ça passe

```


```{r}

M1 <- lmer(Latarena ~ Generation + Replicate + Length + Generation:Replicate + Generation:Length + (1|ID), data=test1, REML=T)
summary(M1)
anova(M1)
ranova(M1) # it's to check the random effect with the likelihood ratio test ; if the random variable is significant, it means that the model is worse without the random effect > keep the random effect.
#On est à O.06 mais c'est mieux de le garder


#Autre méthode : Anova avec un grand "A" -> compare toutes les combinaisons possibles
m1 <- lmer(Latarena ~ Generation*Replicate*Length + (1|ID), data=test1, REML=T)
Anova(m1)

# now we compare the fixed effects with ML, REML = F
M1 <- lmer(Latarena ~ Generation + Replicate + Length + Generation:Replicate + Generation:Length + (1|ID), data=test1, REML=F)
Anova(M1)
M2 <- lmer(Latarena ~ Generation + Replicate + Length + Generation:Length + (1|ID), data=test1, REML=F)
Anova(M2)
M3 <- lmer(Latarena ~ Generation + Replicate + Generation:Length + (1|ID), data=test1, REML=F)
Anova(M3)
M4 <- lmer(Latarena ~ Generation + Replicate + (1|ID), data=test1, REML=F)
Anova(M4)
M5 <- lmer(Latarena ~ Generation + (1|ID), data=test1, REML=F)


anova(M1,M2,M3,M4,M5) # --> M5, pas de diff signif
anova(M5)
ranova(M5)

summary(M5)

```

################################################################################
#                          MODELE VALIDE POUR Latarena : M5
################################################################################

@FS: fort effet generation positif

```{r}

M5 <- lmer(Latarena ~ Generation + (1|ID), data=test1, REML=T)

#Model interpretation
tab <- tab_model(M5, p.val="kr", show.df=T, show.reflvl=T, p.style="scientific_stars")
confint(M5,level=0.95)
tab

#ICC = 0,23


df <- ddply(EPP, c("Generation", "Replicate"), summarize, Mean = mean(Latarena), SD = sd(Latarena))

ggplot(df) +
  aes(x = Generation, y = Mean, fill=Replicate) + 
  scale_fill_manual(breaks = c("1", "2"), values=c("darkgrey","black"))+
  geom_bar(position=position_dodge(0.65), colour="black", stat = "identity", width=0.55)  +
  geom_signif(comparisons = list(c("0","5")), y_position = 1150, tip_length = 1.5, map_signif_level = TRUE, annotation = "***") +
  theme_bw() +
  ylab("Latency to first entry in arena (s)") +
  theme(aspect.ratio = 1/1)+
  scale_y_continuous(limits=c(0,1200), breaks=seq(0, 1300, 200))+
  geom_errorbar(position=position_dodge(0.65), aes(ymin =Mean , ymax = Mean + SD), width = 0.55)

```


```{r}
## check for homogeneity of variance : 
plot(resid(M5)~fitted(M5), xlab="Predicted values", ylab="Normalized residuals")+
  abline(h=0, lty=2)

op <- par(mfrow=c(2,2), mar=c(4,4,.5,.5))

plot(resid(M5)~test1$Generation, xlab="Generation", ylab="Normalized residuals")+
  abline(h=0, lty=2)

plot(resid(M5)~test1$Replicate, xlab="Replicate", ylab="Normalized residuals")+
  abline(h=0, lty=2)

plot(resid(M5)~test1$Length, xlab="Length", ylab="Normalized residuals")+
  abline(h=0, lty=2)

# check for normality of residuals
hist(resid(M5))

qqnorm(resid(M5))
qqline(resid(M5))

```

validation n'est pas parfaite mais c'est ok



Calcul de la répétabilité conditionnelle avec rpt

```{r}

test2 <- test1[test1$Generation=="5",] 
summary(test2)

rpt5 <- rpt(Latarena ~ (1|ID),grname="ID", data= test2, datatype="Gaussian", nboot=1000, npermut=0)
print(rpt5)
#R = 0,244 p=0,112

test3 <- test1[test1$Generation=="0",] 
summary(test3)

rpt0 <- rpt(Latarena ~ (1|ID), grname="ID", data= test3, datatype="Gaussian", nboot=1000, npermut=0)
print(rpt0)

#R=0,217 p=0,0896

#variabilité entre individus plus importante que au sein des ind par rapport aux F5 

#individualité pour F0 mais pas pour F5 (ils perdent leur personnalité avec domestication ?)

?rpt
```

ICC rpz la répétabilité qui est le rapport entre la variabilité entre ind avec la variabilité totale(=variab entre ind et au sein des ind, donc entre 2 réplicats ici)

si les poissons sont + différents entre eux qu'au sein d'un mm poisson ça veut dire qu'on a un ICC élevé et donc qu'ils expriment une personnalité 

au + les poissons sont diff entre eux au + ils expriment une personnalité 


# 2. ----------------- Total time in arena -----------------------------------------------

```{r}

test4 <- EPP[,c("Generation", "ID", "Replicate", "Length", "TTA")] #subset only the columns of interest
test4 <- na.omit(test4) #delete the rows with NA in all the columns

summary(test4)

```

```{r}
test4
# search for extreme values
dotchart(test4$TTA, main="Total time in arena", group=test4$Generation)
dev.off()
hist(test4$TTA)
#visualiser 
#Pas normal, homogene pas suivie
#Pour remédier à ça : on transforme nos données

bestNormalize(test4$TTA) # > orderNorm > prendre Yeo-J
bestNormalize(test4$Length) # > orderNorm > prendre Yeo-J

LA <- yeojohnson(test4$TTA)
x1 <- predict(LA)
test4$TTA <- x1
hist(test4$TTA)

TA <- yeojohnson(test4$Length)
x2 <- predict(TA)
test4$Length <- x2
hist(test4$Length)

```


```{r}
M1 <- lmer(TTA ~ Generation + Replicate + Length + Generation:Replicate + Generation:Length + (1|ID), data=test4, REML=T)
summary(M1)
anova(M1)
ranova(M1) 
#0.059


# now we compare the fixed effects with ML, REML = F
M1 <- lmer(TTA ~ Generation + Replicate + Length + Generation:Replicate + Generation:Length + (1|ID), data=test4, REML=F)
Anova(M1)
M2 <- lmer(TTA ~ Generation + Replicate + Length + Generation:Replicate + (1|ID), data=test4, REML=F)
Anova(M2)
M3 <- lmer(TTA ~ Generation + Replicate + Length + (1|ID), data=test4, REML=F)
Anova(M3)
M4 <- lmer(TTA ~ Generation + Replicate + (1|ID), data=test4, REML=F)
Anova(M4)
M5 <- lmer(TTA ~ Generation + (1|ID), data=test4, REML=F)


anova(M1,M2,M3,M4,M5) #différence quasiment significative pour M4 avec valeur AIC plus faible --> M4
anova(M4)
summary(M4)

M4 <- lmer(TTA ~ Generation + Replicate + (1|ID), data=test4, REML=T)
anova(M4)

tab <- tab_model(M4, p.val="kr", show.df=T, show.reflvl=T, p.style="scientific_stars")
confint(M4,level=0.95)
tab

#fort effet génération et faible effet réplicat


df2 <- ddply(EPP, c("Generation", "Replicate"), summarize, Mean = mean(TTA), SD = sd(TTA))

ggplot(df2) +
  aes(x = Generation, y = Mean, fill=Replicate) + 
  scale_fill_manual(breaks = c("1", "2"), values=c("darkgrey","black"))+
  geom_bar(position=position_dodge(0.65), colour="black", stat = "identity", width=0.55)  +
  geom_signif(comparisons = list(c("0","5")), y_position = 1150, map_signif_level = TRUE, annotation = "***") +
  theme_bw() +
  ylab("Total time in arena (s)") +
  theme(aspect.ratio = 1/1)+
  scale_y_continuous(limits=c(0,1200), breaks=seq(0, 1300, 200))+
  geom_errorbar(position=position_dodge(0.65), aes(ymin =Mean , ymax = Mean + SD), width = 0.55)
   


## check for homogeneity of variance : 
plot(resid(M4)~fitted(M4), xlab="Predicted values", ylab="Normalized residuals")+
  abline(h=0, lty=2)

op <- par(mfrow=c(2,2), mar=c(4,4,.5,.5))

plot(resid(M4)~test4$Generation, xlab="Generation", ylab="Normalized residuals")+
  abline(h=0, lty=2)

plot(resid(M4)~test4$Replicate, xlab="Replicate", ylab="Normalized residuals")+
  abline(h=0, lty=2)

plot(resid(M4)~test4$Length, xlab="Length", ylab="Normalized residuals")+
  abline(h=0, lty=2)

# check for normality of residuals
hist(resid(M4))

qqnorm(resid(M4))
qqline(resid(M4))


test5 <- test4[test4$Generation=="5",] 
summary(test5)

rpt5 <- rpt(TTA ~ Replicate + (1|ID), grname=c("ID", "Fixed"), data= test5, datatype="Gaussian", nboot=1000, npermut=0)
print(rpt5)

#R fixed = 0,012
#R = 0,08 p = 0,336

test6 <- test4[test4$Generation=="0",] 
summary(test6)

rpt0 <- rpt(TTA ~ Replicate + (1|ID), grname=c("ID", "Fixed"), data= test6, datatype="Gaussian", nboot=1000, npermut=0)
print(rpt0)

#R Fixed = 0,028
#R = 0,428 ; p=0,00179


#forte individualité pour F0 mais pas pour F5 (ils perdent leur personnalité avec domestication ?)

```


# 3. ----------------- RTDM -------------------------
  
```{r}

test7 <- EPP[,c("Generation", "ID", "Replicate", "Length", "RTDM")] #subset only the columns of interest
test7 <- na.omit(test7) #delete the rows with NA in all the columns

summary(test7)


# search for extreme values
dotchart(test7$RTDM, main="Relative total distance moved", group=test7$Generation)
dev.off()
hist(test7$RTDM)
#Pas normal, homogeneite pas suivie --> transfo des donnees

bestNormalize(test7$RTDM) # > prendre Yeo-J
bestNormalize(test4$Length) # > Yeo-J

LA <- yeojohnson(test7$RTDM)
x1 <- predict(LA)
test7$RTDM <- x1
hist(test7$RTDM)

TA <- yeojohnson(test7$Length)
x2 <- predict(TA)
test7$Length <- x2
hist(test7$Length)


M1 <- lmer(RTDM ~ Generation + Replicate + Length + Generation:Replicate + Generation:Length + (1|ID), data=test7, REML=T)
summary(M1)
anova(M1)
ranova(M1) 
#0.0218


# now we compare the fixed effects with ML, REML = F
M1 <- lmer(RTDM ~ Generation + Replicate + Length + Generation:Replicate + Generation:Length + (1|ID), data=test7, REML=F)
Anova(M1)
M2 <- lmer(RTDM ~ Generation + Replicate + Length + Generation:Replicate + (1|ID), data=test7, REML=F)
Anova(M2)
M3 <- lmer(RTDM ~ Generation + Replicate + Length + (1|ID), data=test7, REML=F)
Anova(M3)
M4 <- lmer(RTDM ~ Generation + Length + (1|ID), data=test7, REML=F)
Anova(M4)
M5 <- lmer(RTDM ~ Generation + (1|ID), data=test7, REML=F)

anova(M1,M2,M3,M4,M5) # --> M2

M2 <- lmer(RTDM ~ Generation + Replicate + Length + Generation:Replicate + (1|ID), data=test7, REML=T)
anova(M2)

tab <- tab_model(M2, p.val="kr", show.df=T, show.reflvl=T, p.style="scientific_stars")
confint(M2,level=0.95)
tab

#fort effet génération ; fort effet taille (négatif > les petits sont plus actifs) ; l'effet génération dépend du réplicat > habituation différente en fonction de la domestication ; pas d'effet réplicat

df3 <- ddply(test7, c("Generation", "Replicate"), summarize, Mean = mean(RTDM), SD = sd(RTDM))

ggplot(df3) +
  aes(x = Generation, y = Mean, fill=Replicate) + 
  scale_fill_manual(breaks = c("1", "2"), values=c("darkgrey","black"))+
  geom_bar(position=position_dodge(0.65), colour="black", stat = "identity", width=0.55)  +
  geom_signif(comparisons = list(c("0","5")), y_position = 1.3, tip_length = 2.5, vjust = .1, map_signif_level = TRUE, annotation = "***") +
geom_signif(y_position = c(1.2), xmin = c(1.8), 
              xmax = c(2.2), annotation = c("***"),
              tip_length = 0) +
  theme_bw() +
  ylab("Relative total distance moved") +
  theme(aspect.ratio = 1/1)+
  scale_y_continuous(limits=c(0,1.4), breaks=seq(0, 1.4, 0.4))+
  geom_errorbar(position=position_dodge(0.65), aes(ymin =Mean , ymax = Mean + SD), width = 0.55)



ggplot(EPP) +
  aes(x = Replicate, y = RTDM,color=Replicate) + 
  geom_boxplot()+
  theme_bw() +
  facet_grid(~Generation)



plot(EPP$Length, EPP$RTDM)
abline(lm(RTDM~Length, data=EPP),col="red")

ggplot(EPP) +
  aes(x = Length, y = RTDM) + 
  geom_point()+
  theme_bw() +
  geom_smooth()


ggplot(EPP) + 
  aes(x = Length, y = RTDM,color=Replicate)+
  geom_point() +
  theme_bw() +
  ylab("Relative total distance moved") +
  facet_grid(~Generation)+
  geom_smooth()

## check for homogeneity of variance : 
plot(resid(M2)~fitted(M2), xlab="Predicted values", ylab="Normalized residuals")+
  abline(h=0, lty=2)

op <- par(mfrow=c(2,2), mar=c(4,4,.5,.5))

plot(resid(M2)~test7$Generation, xlab="Generation", ylab="Normalized residuals")+
  abline(h=0, lty=2)

plot(resid(M2)~test7$Replicate, xlab="Replicate", ylab="Normalized residuals")+
  abline(h=0, lty=2)

plot(resid(M2)~test7$Length, xlab="Length", ylab="Normalized residuals")+
  abline(h=0, lty=2)

# check for normality of residuals
hist(resid(M2))

qqnorm(resid(M2))
qqline(resid(M2))


test8 <- test7[test7$Generation=="5",] 
summary(test8)

rpt5 <- rpt(RTDM ~ Replicate + Length + (1|ID), grname=c("ID", "Fixed"), data= test8, datatype="Gaussian", nboot=1000, npermut=0)
print(rpt5)

#R fixed = 0,126 / R = 0,177 p = 0,216

test9 <- test7[test7$Generation=="0",] 
summary(test9)

rpt0 <- rpt(RTDM ~ Replicate + Length + (1|ID), grname=c("ID", "Fixed"), data= test9, datatype="Gaussian", nboot=1000, npermut=0)
print(rpt0)

#R = 0.472 ; p=0.000778
#R fixed = 0,366

#forte individualité pour F0 mais pas pour F5 (ils perdent leur personnalité avec domestication ?)


```


# 4. ----------------- IZ ------------------------------

la thigmotaxie doit être testée avec IZ qui est le temps total dans la zone interne divisé par le temps total dans l'arene ; pour que ce soit indépendant de l'audace du poisson

```{r} 

EPP$IZ <- EPP$TTI/EPP$TTA 

test11 <- EPP[,c("Generation", "ID", "Replicate", "Length","TTI","TTA","IZ")] 
test11 <- na.omit(test11) #delete the rows with NA in all the columns

dotchart(test11$IZ, main="Proportion internal zone", group=test11$Generation)
hist(test11$IZ)
hist(test11$Length-hist(log(test11$IZ))

bestNormalize(test11$Length)
Length <- orderNorm(test11$Length)
x3 <- predict(Length)
test11$Length <- x3
hist(test11$Length)

bestNormalize(test11$IZ)
IZ <- orderNorm(test11$IZ)
x4 <- predict(IZ)
test11$IZ <- x4
hist(test11$IZ)


#model 

M1 <- lmer(IZ ~ Generation + Replicate + Length + Generation:Replicate + Generation:Length + (1|ID), data=test11, REML=F)
Anova(M1)

M2 <- lmer(IZ ~ Generation + Replicate + Length + Generation:Replicate + (1|ID), data=test11, REML=F)
Anova(M2)

M3 <- lmer(IZ ~ Generation + Replicate + Length + (1|ID), data=test11, REML=F)
Anova(M3)

M4 <- lmer(IZ ~ Generation + Length + (1|ID), data=test11, REML=F)
Anova(M4)

M5 <- lmer(IZ ~ Generation + (1|ID), data=test11, REML=F)

anova(M1,M2,M3,M4,M5) #--> M5


M5 <- lmer(IZ ~ Generation + (1|ID), data=test11, REML=T)
summary(M5)
anova(M5)


ggplot(EPP) +
  aes(x = Replicate, y = IZ,color=Replicate) + 
  geom_boxplot()+
  theme_bw() +
  facet_grid(~Generation)


#fort effet génération

tab <- tab_model(M5, p.val="kr", show.df=T, show.reflvl=T, p.style="scientific_stars")
confint(M5,level=0.95)
tab



df5 <- ddply(test11, c("Generation", "Replicate"), summarize, Mean = mean(IZ), SD = sd(IZ))

ggplot(df5) +
  aes(x = Generation, y = Mean, fill=Replicate) + 
  scale_fill_manual(breaks = c("1", "2"), values=c("darkgrey","black"))+
  geom_bar(position=position_dodge(0.65), colour="black", stat = "identity", width=0.55)  +
  geom_signif(comparisons = list(c("0","5")), y_position = 0.4, tip_length = 2.5, vjust = .1, map_signif_level = TRUE, annotation = "***") +
  theme_bw() +
  ylab("Relative total time in inner zone") +
  theme(aspect.ratio = 1/1)+
  scale_y_continuous(limits=c(0,0.5), breaks=seq(0, 0.5, 0.1))+
  geom_errorbar(position=position_dodge(0.65), aes(ymin =Mean , ymax = Mean + SD), width = 0.55)



## check for homogeneity of variance : 
plot(resid(M5)~fitted(M5), xlab="Predicted values", ylab="Normalized residuals")+
  abline(h=0, lty=2)

op <- par(mfrow=c(2,2), mar=c(4,4,.5,.5))

plot(resid(M5)~test11$Generation, xlab="Generation", ylab="Normalized residuals")+
  abline(h=0, lty=2)

plot(resid(M5)~test11$Replicate, xlab="Replicate", ylab="Normalized residuals")+
  abline(h=0, lty=2)

plot(resid(M5)~test11$Length, xlab="Length", ylab="Normalized residuals")+
  abline(h=0, lty=2)

# check for normality of residuals
hist(resid(M5))

qqnorm(resid(M5))
qqline(resid(M5))

```

```{r} 

#Repeatability 

M5 <- lmer(IZ ~ Generation + (1|ID), data=test11, REML=T)

test14 <- test11[test11$Generation=="5",] 
summary(test14)

rpt5 <- rpt(IZ ~ (1|ID), grname="ID", data= test14, nboot=1000, npermut=0)
print(rpt5)
#R=0,136 p = 0,275

test15 <- test11[test11$Generation=="0",]
summary(test15)

rpt0 <- rpt(IZ ~ (1|ID), grname="ID", data= test15, nboot=1000, npermut=0)
print(rpt0)

#R=0,45 p=0,00132

#plus d'individualité pour F0 encore 

```