README.md

robotod_stat

packages needed

install.packages('tidyverse', 'glue', 'glue','lme4','multcomp','parameters', 'effectsize','performance','ggpubr')

packages loading

#General utility packages
library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.1.1     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(glue)
library(knitr)
library(data.table)

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

## The following object is masked from 'package:purrr':
## 
##     transpose

#Stats packages
library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(multcomp)

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

library(parameters)

## Registered S3 method overwritten by 'parameters':
##   method                         from      
##   format.parameters_distribution datawizard

library(effectsize)
library(performance)
library(emmeans)
#Plot packages
library(ggpubr)

session information

#output session info
sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## system code page: 65001
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggpubr_0.4.0       emmeans_1.8.1-1    performance_0.8.0  effectsize_0.6.0.1
##  [5] parameters_0.17.0  multcomp_1.4-18    TH.data_1.1-0      MASS_7.3-54       
##  [9] survival_3.2-13    mvtnorm_1.1-3      lme4_1.1-27.1      Matrix_1.3-4      
## [13] data.table_1.14.2  knitr_1.37         glue_1.6.0         forcats_0.5.1     
## [17] stringr_1.4.0      dplyr_1.0.7        purrr_0.3.4        readr_2.1.1       
## [21] tidyr_1.1.4        tibble_3.1.6       ggplot2_3.3.5      tidyverse_1.3.1   
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-153       fs_1.5.2           lubridate_1.8.0    insight_0.17.0    
##  [5] httr_1.4.2         tools_4.1.2        backports_1.4.1    utf8_1.2.2        
##  [9] R6_2.5.1           DBI_1.1.2          colorspace_2.0-2   withr_2.5.0       
## [13] tidyselect_1.1.1   compiler_4.1.2     cli_3.1.1          rvest_1.0.2       
## [17] xml2_1.3.3         sandwich_3.0-1     bayestestR_0.11.5  scales_1.1.1      
## [21] digest_0.6.29      minqa_1.2.4        rmarkdown_2.11     pkgconfig_2.0.3   
## [25] htmltools_0.5.2    dbplyr_2.1.1       fastmap_1.1.0      rlang_0.4.12      
## [29] readxl_1.3.1       rstudioapi_0.13    generics_0.1.1     zoo_1.8-9         
## [33] jsonlite_1.7.3     car_3.0-12         magrittr_2.0.1     Rcpp_1.0.8        
## [37] munsell_0.5.0      fansi_1.0.2        abind_1.4-5        lifecycle_1.0.1   
## [41] stringi_1.7.6      yaml_2.2.1         carData_3.0-5      grid_4.1.2        
## [45] crayon_1.4.2       lattice_0.20-45    haven_2.4.3        splines_4.1.2     
## [49] hms_1.1.1          pillar_1.6.4       boot_1.3-28        estimability_1.4.1
## [53] ggsignif_0.6.3     codetools_0.2-18   reprex_2.0.1       evaluate_0.14     
## [57] modelr_0.1.8       vctrs_0.3.8        nloptr_1.2.2.3     tzdb_0.2.0        
## [61] cellranger_1.1.0   gtable_0.3.0       assertthat_0.2.1   datawizard_0.4.0  
## [65] xfun_0.29          xtable_1.8-4       broom_0.7.11       rstatix_0.7.0     
## [69] ellipsis_0.3.2

loading data

df <- read_csv('https://gitlab.socsci.ru.nl/preclinical-neuroimaging/robotmod/-/raw/main/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(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))

## Rows: 276 Columns: 15

## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr  (2): sex, group
## dbl (13): rat_ID, baseline.day, robot.day, trials, age, restriction, food.di...

## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

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

summary of the data

summary(df) %>% kable("pipe")

	rat_ID	baseline.day	robot.day	trials	sex	group	age	restriction	food.distance.cm	got.food	leaving.sec	total_approaching.sec	backing.sec	foraging.sec	approaching.sec	leaving.score	approaching.score	foraging.score
	407700 : 8	Min. :0.0000	Min. :0.0000	1:60	F:138	LE:156	70 :216	45: 60	25.4 :60	Min. :0.0000	Min. : 1.00	Min. : 2.00	Min. : 1.0	Min. : 4.0	Min. : 0.00	Min. :0.000	Min. :0.000	Min. :0.000
	407701 : 8	1st Qu.:0.0000	1st Qu.:1.0000	2:60	M:138	SD:120	140: 60	60: 96	50.8 :60	1st Qu.:1.0000	1st Qu.: 4.00	1st Qu.: 9.00	1st Qu.: 3.0	1st Qu.: 14.0	1st Qu.: 2.00	1st Qu.:0.000	1st Qu.:0.000	1st Qu.:2.000
	407702 : 8	Median :0.0000	Median :1.0000	3:60				90:120	76.2 :60	Median :1.0000	Median : 7.00	Median : 20.00	Median : 13.0	Median : 35.0	Median : 9.00	Median :1.000	Median :2.000	Median :4.000
	407703 : 8	Mean :0.7826	Mean :0.8696	4:48					101.6:48	Mean :0.7536	Mean :13.21	Mean : 35.01	Mean : 32.8	Mean : 54.2	Mean : 23.67	Mean :2.442	Mean :3.029	Mean :3.935
	407704 : 8	3rd Qu.:0.0000	3rd Qu.:1.0000	5:48					127 :48	3rd Qu.:1.0000	3rd Qu.:13.00	3rd Qu.: 41.00	3rd Qu.: 50.0	3rd Qu.: 80.5	3rd Qu.: 30.00	3rd Qu.:4.000	3rd Qu.:7.000	3rd Qu.:5.250
	407705 : 8	Max. :6.0000	Max. :1.0000							Max. :1.0000	Max. :89.00	Max. :191.00	Max. :262.0	Max. :282.0	Max. :159.00	Max. :7.000	Max. :7.000	Max. :7.000
	(Other):228										NA’s :66	NA’s :78	NA’s :119	NA’s :68	NA’s :78

select testing days (Rday1 here)

test_day <- df %>% filter(df$robot.day == 1)

model comparision

sec/score model comparision (data: SD rat, 90% restrction, foraging, Rday1)

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.

## # Comparison of Model Performance Indices
## 
## Name      |   Model |     AIC | AIC weights |     BIC | BIC weights | R2 (cond.) | R2 (marg.) |   ICC |   RMSE |  Sigma
## -----------------------------------------------------------------------------------------------------------------------
## mod.sec   | lmerMod | 785.240 |     < 0.001 | 813.674 |     < 0.001 |      0.433 |      0.068 | 0.392 | 36.896 | 44.019
## mod.score | lmerMod | 484.861 |        1.00 | 518.311 |        1.00 |      0.634 |      0.064 | 0.609 |  1.227 |  1.413
## [1] "checking model with seconds"
## Warning: Non-normality of residuals detected (p = 0.002).
## Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).
## [1] "checking model with scores"
## OK: residuals appear as normally distributed (p = 0.455).
## OK: Error variance appears to be homoscedastic (p = 0.702).

age effect (data: SD rat, 90% restrction, backing, Rday1)

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)
eta_squared(SD_foraging)

## # Effect Size for ANOVA (Type III)
## 
## Parameter  | Eta2 (partial) |       95% CI
## ------------------------------------------
## sex        |           0.03 | [0.00, 1.00]
## trials     |           0.16 | [0.04, 1.00]
## age        |           0.54 | [0.28, 1.00]
## sex:trials |           0.01 | [0.00, 1.00]
## sex:age    |           0.10 | [0.00, 1.00]
## trials:age |       2.37e-03 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

age effect plot, fig 1a 1b

SD_Rday1_sum<-df.SD_Rday1 %>% group_by(trials, age) %>% summarize(percentage = sum(foraging.score != 7)/n())

## `summarise()` has grouped output by 'trials'. You can override using the `.groups` argument.

Fig1a <-ggbarplot(SD_Rday1_sum, x="trials", y="percentage",
  color = "age", palette = c("#66C2A5", "#FC8D62"),
  size = 2, ylab = "foraging success [%]",xlab = "positions", 
  position = position_dodge(1))

Fig1b<- ggplot(df.SD_Rday1, aes(x=trials, y=foraging.score, group=rat_ID)) + geom_line(color="gray", alpha=0.2)+
stat_summary(fun=mean, aes(group = age, color= age), geom="line", linewidth=2)+
stat_summary(fun.data = mean_se, aes(group=age, color = age), geom = "errorbar", alpha=1, width=0.05)+
scale_color_manual(values = c("#66C2A5", "#FC8D62"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "foraging score") + theme_classic() + facet_wrap(~age, ncol=2)

Fig1a_1b <- ggarrange(Fig1a, Fig1b, ncols=2,labels=c('A','B'), common.legend=TRUE,legend="bottom", widths=c(1,2))
Fig1a_1b

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig1a_1b.svg', plot = Fig1a_1b, device = 'svg',dpi = 300, width = 116, height = 140, units = "mm")

restriciton effect, (data: LE rat, 45%/60% restrction, foraging, Rday1)

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)
eta_squared(LE_foraging)

## # Effect Size for ANOVA (Type III)
## 
## Parameter          | Eta2 (partial) |       95% CI
## --------------------------------------------------
## sex                |           0.33 | [0.07, 1.00]
## trials             |           0.50 | [0.37, 1.00]
## restriction        |           0.33 | [0.07, 1.00]
## sex:trials         |           0.02 | [0.00, 1.00]
## sex:restriction    |           0.11 | [0.00, 1.00]
## trials:restriction |           0.07 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

restriction effect plot, fig 2a 2b

LE_Rday1_sum<-df.LE_Rday1 %>% group_by(trials, restriction) %>% summarize(percentage = sum(foraging.score != 7)/n())

## `summarise()` has grouped output by 'trials'. You can override using the `.groups` argument.

Fig2a <-ggbarplot(LE_Rday1_sum, x="trials", y="percentage",
  color = "restriction", palette = c("#8DA0CB", "#E78AC3"),
  size = 2, ylab = "foraging success [%]",xlab = "positions", 
  position = position_dodge(1))

Fig2b<- ggplot(df.LE_Rday1, aes(x=trials, y=foraging.score, group=rat_ID)) + geom_line(color="gray", alpha=0.2)+
stat_summary(fun=mean, aes(group = restriction, color= restriction), geom="line", linewidth=2)+
stat_summary(fun.data = mean_se, aes(group=restriction, color = restriction), geom = "errorbar", alpha=1, width=0.05)+
scale_color_manual(values = c("#8DA0CB", "#E78AC3"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "foraging score") + theme_classic() + facet_wrap(~restriction, ncol=2,labeller = labeller(group = c("45" = "45%", "60" = "60%")))

Fig2a_2b <- ggarrange(Fig2a, Fig2b, ncols=2,labels=c('A','B'), common.legend=TRUE,legend="bottom", widths=c(1,2))
Fig2a_2b

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig2a_2b.svg', plot = Fig2a_2b, device = 'svg',dpi = 300, width = 116, height = 140, units = "mm")

sex effect, (data: LE and SD rats, 45% and 60% restriction, foraging, Rday1 and Bline day6)

SD rats, PND70

PND70<-df.SD_Rday1 %>% filter(df.SD_Rday1$age == "70")
EFsex_age_70 <- ggplot(PND70, aes(x=trials, y=foraging.score, group=rat_ID))
Fig_SD_sex_age_70 <- EFsex_age_70 +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = sex, color= sex), geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, geom = "errorbar", aes(group=sex, color = sex), alpha=1, width=0.05)+
scale_color_manual(values = c("#66C2A5", "#66C2A5"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "foraging score") + theme_classic() + facet_wrap(~ sex, ncol = 2)
Fig_SD_sex_age_70

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig_SD_sex_age_70.svg', plot = Fig_SD_sex_age_70, device = 'svg',dpi = 300)

## Saving 7 x 5 in image

SD rats, PND140

PND140<-df.SD_Rday1 %>% filter(df.SD_Rday1$age == "140")
EFsex_age_140 <- ggplot(PND140, aes(x=trials, y=foraging.score, group=rat_ID))
Fig_SD_sex_age_140 <- EFsex_age_140 +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1,aes(group = sex, color= sex),geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se,aes(group = sex, color= sex),geom = "errorbar", alpha=1, width=0.05)+
scale_color_manual(values = c("#FC8D62", "#FC8D62"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "foraging score") + theme_classic() + facet_wrap(~ sex, ncol = 2)
Fig_SD_sex_age_140

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig_SD_sex_age_140.svg', plot = Fig_SD_sex_age_140, device = 'svg',dpi = 300)

## Saving 7 x 5 in image

LE rats, 45% restriction

restriction45<-df.LE_Rday1 %>% filter(df.LE_Rday1$restriction == "45")
EFsex_restriction_45 <- ggplot(restriction45, aes(x=trials, y=foraging.score, group=rat_ID))
Fig_LE_sex_restriction_45 <- EFsex_restriction_45 +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = sex, color= sex),geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, aes(group = sex, color= sex),geom = "errorbar", alpha=1, width=0.05)+
scale_color_manual(values = c("#8DA0CB", "#8DA0CB"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "foraging score") + theme_classic() + facet_wrap(~ sex, ncol = 2)
Fig_LE_sex_restriction_45

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig_LE_sex_restriction_45.svg', plot = Fig_LE_sex_restriction_45, device = 'svg',dpi = 300)

## Saving 7 x 5 in image

LE rats, 60% restriction

restriction60<-df.LE_Rday1 %>% filter(df.LE_Rday1$restriction == "60")
EFsex_restriction_60 <- ggplot(restriction60, aes(x=trials, y=foraging.score, group=rat_ID))
Fig_LE_sex_restriction_60 <- EFsex_restriction_60 +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = sex, color= sex), geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, aes(group = sex, color= sex), geom = "errorbar", alpha=1, width=0.05)+
scale_color_manual(values = c("#E78AC3", "#E78AC3"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "foraging score") + theme_classic() + facet_wrap(~ sex, ncol = 2)
Fig_LE_sex_restriction_60

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig_LE_sex_restriction_60.svg', plot = Fig_LE_sex_restriction_60, device = 'svg',dpi = 300)

## Saving 7 x 5 in image

LE rats, pool 45% and 60% restriction data together

p_LE_foraging <- ggplot(df.LE_Rday1, aes(x=trials, y=foraging.score, group=rat_ID))

Fig_LE_sex_general <- p_LE_foraging +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = sex, color= sex), geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, aes(group = sex, color= sex), geom = "errorbar", alpha=1, width=0.05)+
scale_color_manual(values = c("#A6D854", "#FFD92F"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "foraging score") + theme_classic() + facet_wrap(~ sex, ncol = 2)
Fig_LE_sex_general

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig_LE_sex_general.svg', plot = Fig_LE_sex_general, device = 'svg',dpi = 300)

## Saving 7 x 5 in image

LE rats, Blineday6

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=rat_ID))

Fig_LE_sex_blineday6_general <- EFsex_blinedays +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = sex, color= sex), geom = "line", fun.min = 'sd', fun.max = 'sd') +
stat_summary(fun.data = mean_se, aes(group = sex, color= sex), geom = "errorbar", alpha=1, width=0.05)+
scale_color_manual(values = c("#A6D854", "#FFD92F"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "foraging score") + theme_classic() + facet_wrap(~ sex, ncol = 2)

LE_Bday6 <- lmer(foraging.score~ sex + trials + sex:trials + (1|rat_ID),df.LE_Bday6)
eta_squared(LE_Bday6)

Fig_LE_sex_blineday6_general

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig_LE_sex_blineday6_general.svg', plot = Fig_LE_sex_blineday6_general, device = 'svg',dpi = 300)

## Saving 7 x 5 in image

## # Effect Size for ANOVA (Type III)
## 
## Parameter  | Eta2 (partial) |       95% CI
## ------------------------------------------
## sex        |           0.39 | [0.03, 1.00]
## trials     |           0.04 | [0.00, 1.00]
## sex:trials |           0.09 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

sex effect, combine figures, fig 3

Fig3 <- ggarrange(Fig_SD_sex_age_70,Fig_SD_sex_age_140,Fig_LE_sex_restriction_45,Fig_LE_sex_restriction_60,Fig_LE_sex_general,Fig_LE_sex_blineday6_general,ncol=2, nrow=3,labels=c('A','B','C','D','E','F'),common.legend=TRUE,legend="none")
Fig3

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig3.svg', plot = Fig3, device = 'svg',dpi = 300, width = 176, height = 160, units = "mm")

leaving and approaching behaviors (data: LE rats, 45% and 60% restriction, leaving/approaching, Rday1)

mod.score_LE_leaving<- lmer(leaving.score~ restriction + sex + trials + sex:trials + restriction:sex + restriction:trials + (1|rat_ID), df.LE_Rday1)
eta_squared(mod.score_LE_leaving)

mod.score_LE_approaching<- lmer(approaching.score~ restriction + sex + trials +  sex:trials + restriction:sex + restriction:trials + (1|rat_ID), df.LE_Rday1)
eta_squared(mod.score_LE_approaching)

## # Effect Size for ANOVA (Type III)
## 
## Parameter          | Eta2 (partial) |       95% CI
## --------------------------------------------------
## restriction        |           0.04 | [0.00, 1.00]
## sex                |           0.04 | [0.00, 1.00]
## trials             |           0.07 | [0.00, 1.00]
## sex:trials         |           0.05 | [0.00, 1.00]
## restriction:sex    |       6.88e-31 | [0.00, 1.00]
## restriction:trials |           0.02 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).# Effect Size for ANOVA (Type III)
## 
## Parameter          | Eta2 (partial) |       95% CI
## --------------------------------------------------
## restriction        |           0.39 | [0.12, 1.00]
## sex                |           0.02 | [0.00, 1.00]
## trials             |           0.47 | [0.33, 1.00]
## sex:trials         |           0.02 | [0.00, 1.00]
## restriction:sex    |           0.24 | [0.02, 1.00]
## restriction:trials |           0.12 | [0.01, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

here, we further looked into the restriction:sex effect (Eta2 = 0.24 medium size)

mod.score_LE_45_approaching<- lmer(approaching.score~ + sex + trials +  (1|rat_ID), restriction45)
eta_squared(mod.score_LE_45_approaching)

mod.score_LE_60_approaching<- lmer(approaching.score~ + sex + trials +  (1|rat_ID), restriction60)
eta_squared(mod.score_LE_60_approaching)

df.LE_Rday1_approaching_M <- df %>% filter(df$robot.day == 1, df$group == "LE", df$sex == "M")
mod.score_LE_approaching_M<- lmer(approaching.score~ + restriction + trials + (1|rat_ID), df.LE_Rday1_approaching_M)
eta_squared(mod.score_LE_approaching_M)

df.LE_Rday1_approaching_F <- df %>% filter(df$robot.day == 1, df$group == "LE", df$sex == "F")
mod.score_LE_approaching_F<- lmer(approaching.score~ + restriction + trials + (1|rat_ID), df.LE_Rday1_approaching_F)
eta_squared(mod.score_LE_approaching_F)

## # Effect Size for ANOVA (Type III)
## 
## Parameter | Eta2 (partial) |       95% CI
## -----------------------------------------
## sex       |           0.18 | [0.00, 1.00]
## trials    |           0.57 | [0.38, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).# Effect Size for ANOVA (Type III)
## 
## Parameter | Eta2 (partial) |       95% CI
## -----------------------------------------
## sex       |           0.29 | [0.00, 1.00]
## trials    |           0.48 | [0.27, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).# Effect Size for ANOVA (Type III)
## 
## Parameter   | Eta2 (partial) |       95% CI
## -------------------------------------------
## restriction |           0.05 | [0.00, 1.00]
## trials      |           0.44 | [0.23, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).# Effect Size for ANOVA (Type III)
## 
## Parameter   | Eta2 (partial) |       95% CI
## -------------------------------------------
## restriction |           0.64 | [0.26, 1.00]
## trials      |           0.45 | [0.24, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

leaving and approaching plot, fig 4

Fig_LE_leaving <- ggplot(df.LE_Rday1, aes(x=trials, y=leaving.score, group=rat_ID)) + 
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = restriction, color= restriction), geom = "line", fun.min = 'sd', fun.max = 'sd')+
stat_summary(fun.data = mean_se, aes(group = restriction, color= restriction), geom = "errorbar", alpha=1, width=0.05) + 
scale_color_manual(values = c("#8DA0CB", "#E78AC3"))+ 
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "leaving score") + theme_classic() + facet_wrap(~restriction, ncol = 2)
Fig_LE_leaving

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig_LE_leaving.svg', plot = Fig_LE_leaving, device = 'svg',dpi = 300)

## Saving 7 x 5 in image

Fig_LE_approaching <- ggplot(df.LE_Rday1, aes(x=trials, y=approaching.score, group=rat_ID)) +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = restriction, color= restriction), geom = "line", fun.min = 'sd', fun.max = 'sd')+
stat_summary(fun.data = mean_se, aes(group = restriction, color= restriction), geom = "errorbar", alpha=1, width=0.05)+
scale_color_manual(values = c("#8DA0CB", "#E78AC3"))+
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "approaching score")+
theme_classic() + facet_wrap(~ restriction, ncol = 2)
Fig_LE_approaching

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig_LE_approaching.svg', plot = Fig_LE_approaching, device = 'svg',dpi = 300)

## Saving 7 x 5 in image

Fig4 A&B

Fig4 <- ggarrange(Fig_LE_leaving, Fig_LE_approaching, ncols=2,labels=c('A','B'), common.legend=TRUE,legend="none")
Fig4

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig4.svg', plot = Fig4, device = 'svg',dpi = 300, width = 176, height = 140, units = "mm")

Fig4 C&D approaching figure_restriction_sex

fig4c <- ggplot(df.LE_Rday1_approaching_M, aes(x=trials, y=approaching.score, group=rat_ID)) +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = restriction, color= restriction), geom = "line", fun.min = 'sd', fun.max = 'sd')+
stat_summary(fun.data = mean_se, aes(group = restriction, color= restriction), geom = "errorbar", alpha=1, width=0.05) + 
scale_color_manual(values = c("#8DA0CB", "#E78AC3"))+ 
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "approaching score") + theme_classic() + facet_wrap(~restriction, ncol = 2)
fig4c

fig4d <- ggplot(df.LE_Rday1_approaching_F, aes(x=trials, y=approaching.score, group=rat_ID)) +
geom_line(color="gray", alpha=0.2)+
stat_summary(fun = "mean", size = 1, aes(group = restriction, color= restriction), geom = "line", fun.min = 'sd', fun.max = 'sd')+
stat_summary(fun.data = mean_se, aes(group = restriction, color= restriction), geom = "errorbar", alpha=1, width=0.05) + 
scale_color_manual(values = c("#8DA0CB", "#E78AC3"))+ 
scale_y_continuous(breaks = seq(0, 7, by = 1))+
labs(x = "positions", y = "approaching score") + theme_classic() + facet_wrap(~restriction, ncol = 2)
fig4d

Fig4c_d <- ggarrange(Fig_LE_leaving,Fig_LE_approaching,fig4c,fig4d,ncol=2,nrow=2,labels=c('A','B','C','D'),common.legend=TRUE,legend="none")
Fig4c_d

ggsave('C:/Users/xiamen/Desktop/modified robogator/result figs/Fig4c_d.svg', plot = Fig4c_d, device = 'svg',dpi = 300, width = 176, height = 140, units = "mm")