---
title: "robotod_stat"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, 
                      warning = FALSE, 
                      tidy = TRUE, 
                      eval = TRUE, 
                      include = TRUE, 
                      message = TRUE, 
                      results = "hold")
options(knitr.kable.NA = '')
```

## packages needed

```{r, eval=FALSE}
install.packages('tidyverse', 'glue', 'glue','lme4','multcomp','parameters', 'effectsize','performance','ggpubr')
```

## packages loading
```{r}
#General utility packages
library(tidyverse)
library(glue)
library(knitr)
library(data.table)
#Stats packages
library(lme4)
library(multcomp)
library(parameters)
library(effectsize)
library(performance)
library(emmeans)
#Plot packages
library(ggpubr)
```
## session information
```{r sessioninfo, eval=TRUE}
#output session info
sessionInfo()
```

## loading table
```{r}
df <- read_csv('C:/Users/xiamen/Desktop/data_for_R/modified_robogator_paper/robotmod_data.csv') %>% 
          mutate(trials = as.factor(trials)) %>% 
          mutate(group = as.factor(group)) %>% 
          mutate(sex = as.factor(sex)) %>% 
          mutate(rat_ID = as.factor(rat_ID)) %>% 
          mutate(restriction = as.factor(restriction))%>%
          mutate(age = as.factor(age))%>%
          mutate(food.distance.cm = as.factor(food.distance.cm))%>%
         
        # Create score based on time. Following log attribution of score as a function of time. 
          mutate(leaving.score =  case_when(is.na(leaving.sec) ~ 7,
                                                  leaving.sec >= 200 ~ 6,
                                                  leaving.sec >= 100 ~ 5,
                                                  leaving.sec >= 50 ~ 4,
                                                  leaving.sec >= 25  ~ 3,
                                                  leaving.sec >= 12.5 ~ 2,
                                                  leaving.sec >= 6.25 ~ 1,
                                                  leaving.sec <6.25 ~ 0)) %>% 
            mutate(approaching.score =  case_when(is.na(approaching.sec) ~ 7,
                                                  approaching.sec >= 200 ~ 6,
                                                  approaching.sec >= 100 ~ 5,
                                                  approaching.sec >= 50 ~ 4,
                                                  approaching.sec >= 25  ~ 3,
                                                  approaching.sec >= 12.5 ~ 2,
                                                  approaching.sec >= 6.25 ~ 1,
                                                  approaching.sec <6.25 ~ 0)) %>% 
            mutate(backing.score =  case_when(is.na(backing.sec) ~ 7,
                                                  backing.sec >= 200 ~ 6,
                                                  backing.sec >= 100 ~ 5,
                                                  backing.sec >= 50 ~ 4,
                                                  backing.sec >= 25  ~ 3,
                                                  backing.sec >= 12.5 ~ 2,
                                                  backing.sec >= 6.25 ~ 1,
                                                  backing.sec <6.25 ~ 0))  %>%
            mutate(foraging.score =  case_when(is.na(foraging.sec) ~ 7,
                                                  foraging.sec >= 200 ~ 6,
                                                  foraging.sec >= 100 ~ 5,
                                                  foraging.sec >= 50 ~ 4,
                                                  foraging.sec >= 25  ~ 3,
                                                  foraging.sec >= 12.5 ~ 2,
                                                  foraging.sec >= 6.25 ~ 1,
                                                  foraging.sec <6.25 ~ 0))

#head(df) %>% kable("pipe")

```

## summary of the table
```{r}
summary(df) %>% kable("pipe")
```
##select testing days (Rday1 here)
```{r}
test_day <- df %>% filter(df$robot.day == 1)
```

#### model comparision
## sec/score model comparision (data: SD rat, 90% restrction, backing, Rday1) 
```{r}
df.SD_Rday1 <- df %>% filter(df$robot.day == 1, df$group == "SD")
p_sec <- ggplot(df.SD_Rday1, aes(x=trials, y=foraging.sec, group=group, colour=group, shape=sex))
p_sec + stat_summary(fun = "mean", colour = "darkgrey", size = 0.5, geom = "line", fun.min = 'sd', fun.max = 'sd')+
    stat_summary(fun.data = mean_se, geom = "ribbon", alpha=0.3, aes(colour = group, fill = group))+ geom_point(alpha = 1)+theme_minimal()

p_sco <- ggplot(df.SD_Rday1, aes(x=trials, y=foraging.score, group=group, colour=group, shape=sex))
p_sco + stat_summary(fun = "mean", colour = "darkgrey", size = 0.5, geom = "line", fun.min = 'sd', fun.max = 'sd')+
    stat_summary(fun.data = mean_se, geom = "ribbon", alpha=0.3, aes(colour = group, fill = group))+ geom_point(alpha = 1)+theme_minimal()
#***here we noticed,in sec dataset, many missing values in position5, so score dataset is more accurate, besides, the score dataset seems to recapitulate sec dataset

mod.sec <- lmer(foraging.sec~ trials:sex  + trials + sex + (1|rat_ID), df.SD_Rday1)
mod.score <- lmer(foraging.score~ trials:sex  + trials + sex + (1|rat_ID), df.SD_Rday1)
compare_performance(mod.sec, mod.score)
#***The model with scores has a lower AIC and BIC, but comparable R2. This suggest the model with scores is preferable.

print('checking model with seconds')
check_normality(mod.sec)
check_heteroscedasticity(mod.sec)
print('checking model with scores')
check_normality(mod.score)
check_heteroscedasticity(mod.score)
#***Model with scores seem to better follow normal distribution.
```

####age effect, SD rat (70/100 days, M/F, trials 1-5)
##ANOVA-based significant test
```{r}
df.SD_Rday1 <- df %>% filter(df$robot.day == 1, df$group == "SD")
SD_foraging <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID), df.SD_Rday1)

SD_foraging_sex_trials <- update(SD_foraging, . ~ . - sex:trials)
anova(SD_foraging, SD_foraging_sex_trials)

SD_foraging_sex_age <- update(SD_foraging, . ~ . - sex:age)
anova(SD_foraging, SD_foraging_sex_age)

SD_foraging_trials_age <- update(SD_foraging, . ~ . - trials:age)
anova(SD_foraging, SD_foraging_trials_age)

SD_foraging_sex <- lmer(foraging.score~ trials + age + trials:age + (1|rat_ID), df.SD_Rday1)
anova(SD_foraging, SD_foraging_sex)

SD_foraging_trials <- lmer(foraging.score~ sex + age + sex:age + (1|rat_ID), df.SD_Rday1)
anova(SD_foraging, SD_foraging_trials)

SD_foraging_age <- lmer(foraging.score~ sex + trials + sex:trials + (1|rat_ID), df.SD_Rday1)
anova(SD_foraging, SD_foraging_age)

data1_2 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 1|df.SD_Rday1$trials== 2)
data1_3 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 1|df.SD_Rday1$trials== 3)
data1_4 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 1|df.SD_Rday1$trials== 4)
data1_5 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 1|df.SD_Rday1$trials== 5)
data2_3 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 2|df.SD_Rday1$trials== 3)
data2_4 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 2|df.SD_Rday1$trials== 4)
data2_5 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 2|df.SD_Rday1$trials== 5)
data3_4 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 3|df.SD_Rday1$trials== 4)
data3_5 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 3|df.SD_Rday1$trials== 5)
data4_5 <- df.SD_Rday1 %>% filter(df.SD_Rday1$trials == 4|df.SD_Rday1$trials== 5)

SD1_2 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data1_2 )
SD1_3 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data1_3 )
SD1_4 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data1_4 )
SD1_5 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data1_5 )
SD2_3 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data2_3 )
SD2_4 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data2_4 )
SD2_5 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data2_5 )
SD3_4 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data3_4 )
SD3_5 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data3_5)
SD4_5 <- lmer(foraging.score~ sex + trials + age + sex:trials + sex:age + trials:age + (1|rat_ID),data4_5)
print("SD1_2")
anova(SD1_2)
print("SD1_3")
anova(SD1_3)
print("SD1_4")
anova(SD1_4)
print("SD1_5")
anova(SD1_5)
print("SD2_3")
anova(SD2_3)
print("SD2_4")
anova(SD2_4)
print("SD2_5")
anova(SD2_5)
print("SD3_4")
anova(SD3_4)
print("SD3_5")
anova(SD3_5)
print("SD4_5")
anova(SD4_5)

emm_age_trials <- emmeans(SD_foraging, pairwise ~ age | trials)
emm_age_trials
emm_trials_age <- emmeans(SD_foraging, pairwise ~ trials | age)
emm_trials_age
```

##effect size
```{r}
eta_squared(SD_foraging)
```
##effect size for different ages
```{r}
PND70<-df.SD_Rday1 %>% filter(df.SD_Rday1$age == "70")
ESPND70 <- lmer(foraging.score~ sex + trials + sex:trials + (1|rat_ID), PND70)
PND140<-df.SD_Rday1 %>% filter(df.SD_Rday1$age == "140")
ESPND140 <- lmer(foraging.score~ sex + trials + sex:trials + (1|rat_ID), PND140)
eta_squared(ESPND70)
eta_squared(ESPND140)
ESPND70_sex <- lmer(foraging.score~  trials + (1|rat_ID), PND70)
ESPND140_sex <- lmer(foraging.score~ trials  + (1|rat_ID), PND140)
anova(ESPND70,ESPND70_sex)
anova(ESPND140,ESPND140_sex)
```

##plot
```{r}
p_SD_foraging <- ggplot(df.SD_Rday1, aes(x=trials, y=foraging.score, group=age, colour=age))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
p_SD_foraging +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = age, fill = age))+
scale_color_manual(values = colors)+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
#facet_wrap(~ age, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```

## SD rat, success rate
```{r}
#PND70
#males
df.SD_Rday1_70_M <- df.SD_Rday1 %>% filter(df.SD_Rday1$sex == "M",df.SD_Rday1$age == "70")
toselect_70_M <- c("group","food.distance.cm","got.food")
success_rate_70_M <- df.SD_Rday1_70_M %>% dplyr::select(toselect_70_M) %>% data.table() %>% data.table::dcast(food.distance.cm ~ ., fun=list(mean), value.var = "got.food",)
coln<-c("food.distance.cm","PND70.male.success.rate(%)")
colnames(success_rate_70_M) <-coln
#females
df.SD_Rday1_70_F <- df.SD_Rday1 %>% filter(df.SD_Rday1$sex == "F",df.SD_Rday1$age == "70")
toselect_70_F <- c("group","food.distance.cm","got.food")
success_rate_70_F <- df.SD_Rday1_70_F %>% dplyr::select(toselect_70_F) %>% data.table() %>% data.table::dcast(food.distance.cm ~ ., fun=list(mean), value.var = "got.food",)
coln<-c("food.distance.cm","PND70.female.success.rate(%)")
colnames(success_rate_70_F) <-coln
success_rate_70<-left_join(success_rate_70_M,success_rate_70_F)
#PND140
#males
df.SD_Rday1_140_M <- df.SD_Rday1 %>% filter(df.SD_Rday1$sex == "M",df.SD_Rday1$age == "140")
toselect_140_M <- c("group","food.distance.cm","got.food")
success_rate_140_M <- df.SD_Rday1_140_M %>% dplyr::select(toselect_140_M) %>% data.table() %>% data.table::dcast(food.distance.cm ~ ., fun=list(mean), value.var = "got.food",)
coln<-c("food.distance.cm","PND100.male.success.rate(%)")
colnames(success_rate_140_M) <-coln
#females
df.SD_Rday1_140_F <- df.SD_Rday1 %>% filter(df.SD_Rday1$sex == "F",df.SD_Rday1$age == "140")
toselect_140_F <- c("group","food.distance.cm","got.food")
success_rate_140_F <- df.SD_Rday1_140_F %>% dplyr::select(toselect_140_F) %>% data.table() %>% data.table::dcast(food.distance.cm ~ ., fun=list(mean), value.var = "got.food",)
coln<-c("food.distance.cm","PND140.female.success.rate(%)")
colnames(success_rate_140_F) <-coln
success_rate_140<-left_join(success_rate_140_M,success_rate_140_F)
success_rate<-left_join(success_rate_70,success_rate_140)
print(success_rate)
```

####restriciton effect, LE rat (45%/60%, M/F, trials 1-5)
## ANOVA-based significant test
```{r}
df.LE_Rday1 <- df %>% filter(df$robot.day == 1, df$group == "LE")

LE_foraging <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID), df.LE_Rday1)

LE_foraging_sex_trials <- update(LE_foraging, . ~ . - sex:trials)
anova(LE_foraging, LE_foraging_sex_trials)

LE_foraging_sex_restriction <- update(LE_foraging, . ~ . - sex:restriction)
anova(LE_foraging, LE_foraging_sex_restriction)

LE_foraging_trials_restriction <- update(LE_foraging, . ~ . - trials:restriction)
anova(LE_foraging, LE_foraging_trials_restriction)

LE_foraging_sex <- lmer(foraging.score~ trials + restriction + trials:restriction + (1|rat_ID), df.LE_Rday1)
anova(LE_foraging, LE_foraging_sex)

LE_foraging_trials <- lmer(foraging.score~ sex + restriction + sex:restriction + (1|rat_ID), df.LE_Rday1)
anova(LE_foraging, LE_foraging_trials)

LE_foraging_restriction <- lmer(foraging.score~ sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1)
anova(LE_foraging, LE_foraging_restriction)

LE_data1_2 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 1|df.LE_Rday1$trials== 2)
LE_data1_3 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 1|df.LE_Rday1$trials== 3)
LE_data1_4 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 1|df.LE_Rday1$trials== 4)
LE_data1_5 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 1|df.LE_Rday1$trials== 5)
LE_data2_3 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 2|df.LE_Rday1$trials== 3)
LE_data2_4 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 2|df.LE_Rday1$trials== 4)
LE_data2_5 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 2|df.LE_Rday1$trials== 5)
LE_data3_4 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 3|df.LE_Rday1$trials== 4)
LE_data3_5 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 3|df.LE_Rday1$trials== 5)
LE_data4_5 <- df.LE_Rday1 %>% filter(df.LE_Rday1$trials == 4|df.LE_Rday1$trials== 5)

LE1_2 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data1_2 )
LE1_3 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data1_3 )
LE1_4 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data1_4 )
LE1_5 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data1_5 )
LE2_3 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data2_3 )
LE2_4 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data2_4 )
LE2_5 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data2_5 )
LE3_4 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data3_4 )
LE3_5 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data3_5)
LE4_5 <- lmer(foraging.score~ sex + trials + restriction + sex:trials + sex:restriction + trials:restriction + (1|rat_ID),LE_data4_5)
print("LE1_2")
anova(LE1_2)
print("LE1_3")
anova(LE1_3)
print("LE1_4")
anova(LE1_4)
print("LE1_5")
anova(LE1_5)
print("LE2_3")
anova(LE2_3)
print("LE2_4")
anova(LE2_4)
print("LE2_5")
anova(LE2_5)
print("LE3_4")
anova(LE3_4)
print("LE3_5")
anova(LE3_5)
print("LE4_5")
anova(LE4_5)

emm_restriction_trials <- emmeans(LE_foraging, pairwise ~ restriction | trials)
emm_restriction_trials
emm_trials_restriction <- emmeans(LE_foraging, pairwise ~ trials | restriction)
emm_trials_restriction

emm_sex_trials <- emmeans(LE_foraging, pairwise ~ sex | trials)
emm_sex_trials
emm_trials_sex <- emmeans(LE_foraging, pairwise ~ trials | sex)
emm_trials_sex
```

##effect size
```{r}
eta_squared(LE_foraging)
```
##effect size for different restrictions
```{r}
restriction45<-df.LE_Rday1 %>% filter(df.LE_Rday1$restriction == "45")
ESrestriction45 <- lmer(foraging.score~ sex + trials + sex:trials + (1|rat_ID), restriction45)
restriction60<-df.LE_Rday1 %>% filter(df.LE_Rday1$restriction == "60")
ESrestriction60 <- lmer(foraging.score~ sex + trials + sex:trials + (1|rat_ID), restriction60)
eta_squared(ESrestriction45)
eta_squared(ESrestriction60)
ESrestriction45_sex <- lmer(foraging.score~  trials + (1|rat_ID), restriction45)
ESrestriction60_sex <- lmer(foraging.score~ trials  + (1|rat_ID), restriction60)
anova(ESrestriction45,ESrestriction45_sex)
anova(ESrestriction60,ESrestriction60_sex)
```

##plot-restriction
```{r}
p_LE_foraging <- ggplot(df.LE_Rday1, aes(x=trials, y=foraging.score, group=restriction, colour=restriction))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
p_LE_foraging +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = restriction, fill = restriction))+
scale_color_manual(values = colors)+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
#facet_wrap(~ age, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```
##plot_sex effect
#PND70
```{r}
EFsex_age_70 <- ggplot(PND70, aes(x=trials, y=foraging.score, group=sex, colour=sex))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
EFsex_age_70 +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = age, fill = age))+
scale_color_manual(values = colors)+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
#facet_wrap(~ age, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```
#PND140
```{r}
EFsex_age_140 <- ggplot(PND140, aes(x=trials, y=foraging.score, group=sex, colour=sex))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
EFsex_age_140 +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = age, fill = age))+
scale_color_manual(values = colors)+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
#facet_wrap(~ age, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```
#45% restriction
```{r}
EFsex_restriction_45 <- ggplot(restriction45, aes(x=trials, y=foraging.score, group=sex, colour=sex))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
EFsex_restriction_45 +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = age, fill = age))+
scale_color_manual(values = colors)+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
#facet_wrap(~ age, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```
#60% restriction
```{r}
EFsex_restriction_60 <- ggplot(restriction60, aes(x=trials, y=foraging.score, group=sex, colour=sex))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
EFsex_restriction_60 +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = age, fill = age))+
scale_color_manual(values = colors)+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
#facet_wrap(~ age, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```
##general sex effect
```{r}
p_LE_backing <- ggplot(df.LE_Rday1, aes(x=trials, y=foraging.score, group=sex, colour=sex))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
p_LE_backing +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = sex, fill = sex))+
scale_color_manual(values = colors)+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
#facet_wrap(~ age, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```
##Blineday6 sex effect plot
```{r}
df.LE_Bday6 <- df %>% filter(df$baseline.day == 6, df$group == "LE")
EFsex_blinedays <- ggplot(df.LE_Bday6, aes(x=trials, y=foraging.score, group=sex, colour=sex))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
EFsex_blinedays +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = sex, fill = sex))+
scale_color_manual(values = colors)+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
#facet_wrap(~ age, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```

## LE rat, success rate
```{r}
#45% food restriction
#males
df.LE_Rday1_45_M <- df.LE_Rday1 %>% filter(df.LE_Rday1$sex == "M",df.LE_Rday1$restriction == "45")
toselect_45_M <- c("group","food.distance.cm","got.food")
success_rate_45_M <- df.LE_Rday1_45_M %>% dplyr::select(toselect_45_M) %>% data.table() %>% data.table::dcast(food.distance.cm ~ ., fun=list(mean), value.var = "got.food",)
coln<-c("food.distance.cm","45%.male.success.rate(%)")
colnames(success_rate_45_M) <-coln
#females
df.LE_Rday1_45_F <- df.LE_Rday1 %>% filter(df.LE_Rday1$sex == "F",df.LE_Rday1$restriction == "45")
toselect_45_F <- c("group","food.distance.cm","got.food")
success_rate_45_F <- df.LE_Rday1_45_F %>% dplyr::select(toselect_45_F) %>% data.table() %>% data.table::dcast(food.distance.cm ~ ., fun=list(mean), value.var = "got.food",)
coln<-c("food.distance.cm","45%.female.success.rate(%)")
colnames(success_rate_45_F) <-coln
success_rate_45<-left_join(success_rate_45_M,success_rate_45_F)
#60% food restriction
#males
df.LE_Rday1_60_M <- df.LE_Rday1 %>% filter(df.LE_Rday1$sex == "M",df.LE_Rday1$restriction == "60")
toselect_60_M <- c("group","food.distance.cm","got.food")
success_rate_60_M <- df.LE_Rday1_60_M %>% dplyr::select(toselect_60_M) %>% data.table() %>% data.table::dcast(food.distance.cm ~ ., fun=list(mean), value.var = "got.food",)
coln<-c("food.distance.cm","60%.male.success.rate(%)")
colnames(success_rate_60_M) <-coln
#females
df.LE_Rday1_60_F <- df.LE_Rday1 %>% filter(df.LE_Rday1$sex == "F",df.LE_Rday1$restriction == "60")
toselect_60_F <- c("group","food.distance.cm","got.food")
success_rate_60_F <- df.LE_Rday1_60_F %>% dplyr::select(toselect_60_F) %>% data.table() %>% data.table::dcast(food.distance.cm ~ ., fun=list(mean), value.var = "got.food",)
coln<-c("food.distance.cm","60%.female.success.rate(%)")
colnames(success_rate_60_F) <-coln
success_rate_60<-left_join(success_rate_60_M,success_rate_60_F)
success_rate<-left_join(success_rate_45,success_rate_60)
print(success_rate)
```

## distance effect, using LE 60% Bline day6
```{r}
df.LE_Bday6 <- df %>% filter(df$baseline.day == 6, df$group == "LE")
LE_Bday6_foraging<- lmer(foraging.score~  sex + food.distance.cm + sex:food.distance.cm + (1|rat_ID), df.LE_Bday6)

LE_Bday6_distance_sex<-lmer(foraging.score~  sex + food.distance.cm + (1|rat_ID), df.LE_Bday6)
LE_Bday6_distance<-lmer(foraging.score~  sex +  (1|rat_ID), df.LE_Bday6)

anova(LE_Bday6_distance_sex, LE_Bday6_foraging)
anova(LE_Bday6_distance, LE_Bday6_foraging)

eta_squared(LE_Bday6_foraging)

```
## distance effect figure
```{r}
p_LE_distance <- ggplot(df.LE_Bday6, aes(x=food.distance.cm, y=foraging.score, group=sex, colour=sex))
colors <- c("#c6c7c9", "#d15034", "#0c0d0f", "#f1ae2d")
p_LE_distance +
stat_summary(fun = "mean", size = 1, geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", alpha=1, width=0.05, aes(colour = sex, fill = sex))+
scale_color_manual(values = colors)+

#facet_wrap(~ food.distance.cm, ncol = 2)
geom_jitter(alpha = 0.5,width = 0.1) + labs(x = "positions", y = "foraging score") + theme_minimal()
```

#------------------------------------
#### some other behavioral parameters (LE rat)
## leaving
```{r}
mod.score_LE_leaving<- lmer(leaving.score~ restriction + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1)
eta_squared(mod.score_LE_leaving)
```
## 45% leaving
```{r}
mod.score_LE_45_leaving<- lmer(leaving.score~ + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1_45)
eta_squared(mod.score_LE_45_leaving)
```
## 60% leaving
```{r}
mod.score_LE_60_leaving<- lmer(leaving.score~ + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1_60)
eta_squared(mod.score_LE_60_leaving)
```
## leaving figure
```{r}
LE_Rday1_leaving <- ggplot(df.LE_Rday1, aes(x=food.distance.cm, y=leaving.score, group=restriction, colour=restriction, shape=sex))
LE_Rday1_leaving + stat_summary(fun = "mean", colour = "darkgrey", size = 0.5, geom = "line", fun.min = 'sd', fun.max = 'sd')+
    stat_summary(fun.data = mean_se, geom = "ribbon", alpha=0.3, aes(colour = restriction, fill = restriction))+
      geom_jitter(alpha = 1,width = 0.1)+theme_minimal()
```

## approaching
```{r}
mod.score_LE_approaching<- lmer(approaching.score~ restriction + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1)
eta_squared(mod.score_LE_approaching)
```
## 45% approaching
```{r}
mod.score_LE_45_approaching<- lmer(approaching.score~ + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1_45)
eta_squared(mod.score_LE_45_approaching)
```
## 60% approaching
```{r}
mod.score_LE_60_approaching<- lmer(approaching.score~ + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1_60)
eta_squared(mod.score_LE_60_approaching)
```
## approaching figure
```{r}
LE_Rday1_approaching <- ggplot(df.LE_Rday1, aes(x=food.distance.cm, y=approaching.score, group=restriction, colour=restriction, shape=sex))
LE_Rday1_approaching + stat_summary(fun = "mean", colour = "darkgrey", size = 0.5, geom = "line", fun.min = 'sd', fun.max = 'sd')+
    stat_summary(fun.data = mean_se, geom = "ribbon", alpha=0.3, aes(colour = restriction, fill = restriction))+
      geom_jitter(alpha = 1,width = 0.1)+theme_minimal()
```
####################
## backing
```{r}
mod.score_LE_backing<- lmer(backing.score~ restriction + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1)
eta_squared(mod.score_LE_backing)
```
## 45% backing
```{r}
mod.score_LE_45_backing<- lmer(backing.score~ + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1_45)
eta_squared(mod.score_LE_45_backing)
```
## 60% backing
```{r}
mod.score_LE_60_backing<- lmer(backing.score~ + sex + trials + sex:trials + (1|rat_ID), df.LE_Rday1_60)
eta_squared(mod.score_LE_60_backing)
```
## backing figure
```{r}
LE_Rday1_backing <- ggplot(df.LE_Rday1, aes(x=food.distance.cm, y=backing.score, group=restriction, colour=restriction, shape=sex))
LE_Rday1_backing + stat_summary(fun = "mean", colour = "darkgrey", size = 0.5, geom = "line", fun.min = 'sd', fun.max = 'sd')+
    stat_summary(fun.data = mean_se, geom = "ribbon", alpha=0.3, aes(colour = restriction, fill = restriction))+
      geom_jitter(alpha = 1,width = 0.1)+theme_minimal()
```
###################