Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
---
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)
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
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_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
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
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)
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
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_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)
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
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()
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
## 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()
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
```
#------------------------------------
#### 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()
```
###################