--- title: "hecale_analysis" author: "Wyatt Toure" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ## R Markdown ```{r} library(lme4) library(lmerTest) library(psych) library(ggplot2) library(ggpubr) library(DHARMa) library(dplyr) library(tidyr) library(DHARMa) library(effects) # REPLACE PATH IN read.csv(PATH) with path to full_heliconius_dataset data = read.csv('PATH') data = data %>% filter(species == 'hecale') data = data %>% group_by(id) %>% mutate(total.activity = (n.correct.training+n.incorrect.training)) ``` ```{r} # Is there a difference in activity between morning and afternoon activity.model = glmer((n.correct.training+n.incorrect.training) ~ time.of.day + (1 | cage/id), data= data, family = poisson) summary(activity.model) data = data %>% group_by(id) %>% mutate(total.activity = (n.correct.training+n.incorrect.training)) data %>% group_by(time.of.day) %>% summarise(mean(total.activity)) ``` ```{r} # Model checking simulationOutput <- simulateResiduals(fittedModel = activity.model) testDispersion(simulationOutput) testUniformity(simulationOutput) testZeroInflation(simulationOutput) # model is zero inflated but we know it is since we are demonstrating that a good chunk of # individuals have low activity in the evening preventing all colour choice combinations from # being made # model ok otherwise ``` ```{r} ### Group by id and remove any individuals who did not make both colour choices in either ### morning or afternoon during training data2 = data %>% group_by(id) %>% filter(all(n.correct.training > 0 & n.incorrect.training > 0)) ### Remove individuals who did not respond at all during a session in final test data3 = data2 %>% group_by(id) %>% filter(all(n.morning.colour.test > 0 | n.afternoon.colour.test > 0)) ### Is there a shift in preference for colours when individuals are naive naive.model = glmer(cbind(n.morning.colour.initial, n.afternoon.colour.initial) ~ time.of.day + (1 | cage/id), data= data3, family = binomial) summary(naive.model) ### No p-value = 0.390 ``` ```{r} ### model checking simulationOutput <- simulateResiduals(fittedModel = naive.model) testDispersion(simulationOutput) testUniformity(simulationOutput = simulationOutput) testZeroInflation(simulationOutput) # model OK ``` ```{r} ### Is there a shift in preference for colours when individuals are trained trained.model = glmer(cbind(n.morning.colour.test, n.afternoon.colour.test) ~ time.of.day * final.presentation + (1 | cage/id), data= data3, family = binomial) summary(trained.model) ### no p = 0.06 for time of day effect ``` ```{r} # model checking simulationOutput <- simulateResiduals(fittedModel = trained.model) testDispersion(simulationOutput) testUniformity(simulationOutput = simulationOutput) testZeroInflation(simulationOutput) # model ok ``` ```{r} ### Remove individuals who did not achieve greater than 50% correct responses ### in last 2 days of training data4 = data3 %>% group_by(id) %>% filter(all(prop.correct.training.last.2 > 0.5)) ### Do individuals who are perform well in last 2 days of training shift preference initially naive.model.high.performers = glmer(cbind(n.morning.colour.initial, n.afternoon.colour.initial) ~ time.of.day + (1 | cage/id), data= data4, family = binomial) summary(naive.model.high.performers) ### no p = 0.79 ``` ```{r} # model checking simulationOutput <- simulateResiduals(fittedModel = naive.model.high.performers) testDispersion(simulationOutput) testUniformity(simulationOutput = simulationOutput) testZeroInflation(simulationOutput) # model ok ``` ```{r} # pull out individuals below criterion data5 = data3 %>% group_by(id) %>% filter(any(prop.correct.training.last.2 < 0.51)) ### Do individuals who were incorrect at high rate shift preference initially ### this model fails to converge, rendering p-values meaningless ### unless we include an observation level random effect this is what ### the next lines of code before the model does df = data.frame(OLRE=c(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)) data5 <- cbind(data5, OLRE = df$OLRE) naive.model.low.performers = glmer(cbind(n.morning.colour.initial, n.afternoon.colour.initial) ~ time.of.day + (1 | cage/id) + (1 | OLRE), data =data5, family = binomial) summary(naive.model.low.performers) ### no p = 0.56 ``` ```{r} # model checking simulationOutput <- simulateResiduals(fittedModel = naive.model.low.performers) testDispersion(simulationOutput) testUniformity(simulationOutput = simulationOutput) testZeroInflation(simulationOutput) # model ok ``` ```{r} ### Do individuals who were correct at a high rate shift preference when trained ### Cage random effect dropped singular fit issue trained.model.high.performers = glmer(cbind(n.morning.colour.test, n.afternoon.colour.test) ~ time.of.day * final.presentation + (1 | id), data= data4, family = binomial) summary(trained.model.high.performers) ### yes p = 0.01 and no significant effect of final presentation order ``` ```{r} # model checking simulationOutput <- simulateResiduals(fittedModel = trained.model.high.performers) testDispersion(simulationOutput) testUniformity(simulationOutput = simulationOutput) testZeroInflation(simulationOutput) # model ok ``` ```{r} ### Do individuals who were not correct at a high rate shift preference when trained trained.model.low.performers = glmer(cbind(n.morning.colour.test, n.afternoon.colour.test) ~ time.of.day * final.presentation + (1 | cage/id), data= data5, family = binomial) summary(trained.model.low.performers) ### no p = 0.24 ``` ```{r} # model checking simulationOutput <- simulateResiduals(fittedModel = trained.model.low.performers) testDispersion(simulationOutput) testUniformity(simulationOutput = simulationOutput) testZeroInflation(simulationOutput) # model ok ``` ```{r} ### Naive preference data plot for high performers g <- ggplot(data4, aes(time.of.day, (n.morning.colour.initial/(n.morning.colour.initial+n.afternoon.colour.initial)))) + theme_classic() a <- g + geom_jitter(size =2, alpha = .25, width = 0.02) + stat_summary( geom = 'point', fun = 'mean', size = 4.5, shape = 15) + stat_summary(geom = 'errorbar', fun.data = 'mean_ci', position = position_dodge(width=0), width = 0.1) naive.pref.plot <- a + ylab('Morning reward colour preference') + xlab("Time of Day") + ggtitle("Naive shift in preference") + theme(axis.text=element_text(size=10), axis.title=element_text(size=11,face="bold"), plot.title = element_text(size=12, hjust=0.5)) + ylim(-0.05, 1.05) + geom_line(aes(group = id), color = 'grey') ### Trained preference data plot for high performers g2 <- ggplot(data4, aes(time.of.day, (n.morning.colour.test/(n.morning.colour.test+n.afternoon.colour.test)))) + theme_classic()# + scale_colour_manual(values = c('#353b48', '#e84118')) b <- g2 + geom_jitter(size =2, alpha = .25, width = 0.02) + stat_summary( geom = 'point', fun = 'mean', size = 4.5, shape = 15) + stat_summary(geom = 'errorbar', fun.data = 'mean_ci', position = position_dodge(width=.05), width = 0.1) trained.pref.plot <- b + ylab('Morning reward colour preference') + xlab("Time of Day") + ggtitle("Trained shift in preference") + theme(axis.text=element_text(size=10), axis.title=element_text(size=11,face="bold"), plot.title = element_text(size=12, hjust=0.5)) + ylim(-0.05, 1.05) + geom_line(aes(group = id), color = 'grey') ggarrange(naive.pref.plot, trained.pref.plot, nrow = 1) #### Plot end ####### ``` Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.