Preliminaries

Load libraries

library(knitr)       #for R Markdown
library(MASS)        #for negative binomial glm (glm.nb)
library(lme4)        #for mixed models
library(emmeans)     #for posthoc
library(car)         #for Anova
library(survival)    #for survival analysis
library(coxme)
library(rptR)        #for repeatability analysis
library(MuMIn)       #for model selection (dredge)
library(ggfortify)   #for plotting survival analysis
library(ggsignif)    #for labeling significance in ggplots
library(GGally)     #for correlation matrix
library(tidyverse)   #for data processing, put last to avoid function masking

Load data and clean up

data_raw <- read.csv("data_maze_040326.csv",                 
                 # to avoid reading errors
                 fileEncoding="UTF-8-BOM", na.strings = "")


data <- data_raw %>%
  
  #drop trailing NA reads
  filter(!if_all(everything(), is.na)) %>%
  
  #correct variable types
  mutate_at(c("ID","sex","assayer","trial","shelter_YN","shelter_quadrant"), as.factor) %>%
  mutate_at(c("found_date","morph_date","assay_date"), lubridate::mdy) %>%
  
  #derive variables
  mutate(
    
    #life history variables
    lat_meta = as.numeric(morph_date - found_date),
    age = as.numeric(assay_date - morph_date),
         
    #performance variables     
    lat_shelter_from_movement = as.numeric(lat_shelter - lat_move),
    lines_rate = lines_crossed/lat_trial,
    grids_rate = grids_explored/lat_trial,
    
    #for survival analysis and plotting
    shelter_10 = recode(shelter_YN, "Y" = 1, "N" = 0),
    move_10 = ifelse(lat_move==lat_trial, 0, 1)
    
    )

Make an aggregrate version of the dataset by frog ID

data_ind <- data %>%
  group_by(ID) %>%
  summarise(lat_move_avg = mean(lat_move),
            lat_shelter_avg = mean(lat_shelter),
            lat_shelter_mov_avg = mean(lat_shelter_from_movement),
            lat_shelter_improv = lat_shelter[trial == 1] - lat_shelter[trial == 5],
            lines_rate_avg = mean(lines_rate),
            grids_rate_avg = mean(grids_rate),
            prop_complete = sum(shelter_YN == "Y")/n(),
            found_date = mean(found_date),
            morph_date = mean(morph_date), 
            age = mean(age),
            lat_meta = mean(lat_meta)
            )

1. Repeatability

  • rpt{rptR}: Repeatability Estimation
  • P_permut: p-value from permutation test (more robust?)
  • LRT_P: p-value from likelihood ratio test

All three personality axes are repeatable!

Boldness

rpt_result <- rpt(lat_move ~ (1 | ID), 
                  grname = "ID", 
                  data = data, 
                  datatype = "Gaussian", 
                  nboot = 1000, 
                  npermut = 1000)
## Bootstrap Progress:
## Permutation Progress for ID :
summary(rpt_result)$rpt %>% kable(digits = 3)
R SE 2.5% 97.5% P_permut LRT_P
rpt 0.218 0.131 0 0.478 0.016 0.032

Activity

rpt_result <- rpt(lines_rate ~ (1 | ID), 
                  grname = "ID", 
                  data = data, 
                  datatype = "Gaussian", 
                  nboot = 1000, 
                  npermut = 1000)
## Bootstrap Progress:
## Permutation Progress for ID :
summary(rpt_result)$rpt %>% kable(digits = 3)
R SE 2.5% 97.5% P_permut LRT_P
rpt 0.451 0.142 0.138 0.673 0.001 0

Exploration

rpt_result <- rpt(grids_rate ~ (1 | ID), 
                  grname = "ID", 
                  data = data, 
                  datatype = "Gaussian", 
                  nboot = 1000, 
                  npermut = 1000)
## Bootstrap Progress:
## Permutation Progress for ID :
summary(rpt_result)$rpt %>% kable(digits = 3)
R SE 2.5% 97.5% P_permut LRT_P
rpt 0.412 0.141 0.099 0.652 0.001 0

2. Behavioral Syndrome

Stats

boldness-activity

mod_boldness_activity <- lmer(lines_rate ~ lat_move + (1|ID),
                              data = data)

summary(mod_boldness_activity)$coefficients
##                  Estimate   Std. Error    t value
## (Intercept)  2.224513e-02 2.225872e-03  9.9938964
## lat_move    -1.150276e-05 1.265838e-05 -0.9087068
Anova(mod_boldness_activity)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: lines_rate
##           Chisq Df Pr(>Chisq)
## lat_move 0.8257  1     0.3635

boldness-exploration

mod_boldness_exploration <- lmer(grids_rate ~ lat_move + (1|ID),
                              data = data)

summary(mod_boldness_exploration)$coefficients
##                 Estimate   Std. Error   t value
## (Intercept) 7.529567e-03 6.981217e-04 10.785465
## lat_move    5.257119e-06 3.647541e-06  1.441277
Anova(mod_boldness_exploration)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: grids_rate
##           Chisq Df Pr(>Chisq)
## lat_move 2.0773  1     0.1495

activity-exploration

mod_activity_exploration <- lmer(grids_rate ~ lines_rate + (1|ID),
                              data = data)

summary(mod_activity_exploration)$coefficients
##                Estimate   Std. Error  t value
## (Intercept) 0.005546651 0.0008851838 6.266101
## lines_rate  0.110757968 0.0354561117 3.123805
Anova(mod_activity_exploration)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: grids_rate
##             Chisq Df Pr(>Chisq)   
## lines_rate 9.7582  1   0.001785 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation plots - by trial

fig_boldness_activity <- 
  ggplot(data, aes(x = lat_move, y = lines_rate)) +
  geom_point(color = "black", size = 2) +
  geom_smooth(method = lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  
  xlab("Latency to move (s)") +
  ylab("# Lines crossed / Trial time (s)") +
  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)

fig_boldness_exploration <- 
  ggplot(data, aes(x = lat_move, y = grids_rate)) +
  geom_point(color = "black", size = 2) +
  geom_smooth(method = lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  
  xlab("Latency to move (s)") +
  ylab("# Grids explored / Trial time (s)") +
  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)

fig_activity_exploration <- 
  ggplot(data, aes(x = lines_rate, y = grids_rate)) +
  geom_point(color = "black", size = 2) +
  geom_smooth(method = lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  
  xlab("# Lines crossed / Trial time (s)") +
  ylab("# Grids explored / Trial time (s)") +
  
  theme_bw(base_size = 15) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)
egg::ggarrange(fig_boldness_activity, fig_boldness_exploration, fig_activity_exploration,
               nrow = 1,
               labels = c("A","B","C"))

Correlation Plots - by individual

fig_boldness_activity_ind <- 
  ggplot(data_ind, aes(x = lat_move_avg, y = lines_rate_avg)) +
  
  # for overlaying trial level data - doesn't look good on the plot though
  # geom_point(data = data, aes(x = lat_move, y = lines_rate),
  #            color = "grey80", size = 1) +
  
  geom_point(color = "black", size = 2) +
  geom_smooth(method = lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  
  xlab("Latency to move (s)") +
  ylab("# Lines crossed / Trial time (s)") +
  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)


fig_boldness_exploration_ind <- 
  ggplot(data_ind, aes(x = lat_move_avg, y = grids_rate_avg)) +
  geom_point(color = "black", size = 2) +
  geom_smooth(method = lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  
  xlab("Latency to move (s)") +
  ylab("# Grids explored / Trial time (s)") +
  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)


fig_activity_exploration_ind <- 
  ggplot(data_ind, aes(x = lines_rate_avg, y = grids_rate_avg)) +
  geom_point(color = "black", size = 2) +
  geom_smooth(method = lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  
  xlab("# Lines crossed / Trial time (s)") +
  ylab("# Grids explored / Trial time (s)") +
  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)
egg::ggarrange(fig_boldness_activity_ind, fig_boldness_exploration_ind, fig_activity_exploration_ind,
               nrow = 1,
               labels = c("A","B","C"))

3. Maze Performance and Life History

stats

mod<- lm(prop_complete ~ lat_meta, data = data_ind)

summary(mod)$coefficients
##                Estimate Std. Error    t value  Pr(>|t|)
## (Intercept) -0.51314286 1.10907807 -0.4626751 0.6535013
## lat_meta     0.01314286 0.01866657  0.7040853 0.4974475
mod<- lm(prop_complete ~ age, data =  data_ind)

summary(mod)$coefficients
##                 Estimate  Std. Error   t value  Pr(>|t|)
## (Intercept) 0.1186611406 0.613392068 0.1934507 0.8504788
## age         0.0003952963 0.001630547 0.2424317 0.8133454

figure

fig_meta_comp <-
  ggplot(data_ind ,
         aes(x = lat_meta, y = prop_complete)) +
  geom_point(color = "black", size = 2) +
  #geom_smooth(method = lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  
  xlab("Latency to Metamorphose (days)") +
  ylab("Percentage of Success (%)") +
  scale_y_continuous(labels = scales::percent) +  
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1)


fig_age_comp <-
  ggplot(data_ind,
         aes(x = age, y = prop_complete)) +
  geom_point(color = "black", size = 2) +

  # geom_smooth(method = lm, alpha = 0.2, color = "SteelBlue", fill = "SteelBlue")+
  
  xlab("Age (days)") +
  ylab("Percentage of Success (%)") +
  scale_y_continuous(labels = scales::percent) +
  theme_bw(base_size = 15)+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        aspect.ratio=1) 
egg::ggarrange(fig_meta_comp, fig_age_comp,
               nrow = 1,
               labels = c("A","B"))

Additional - survival plot

No clear sign of learning! They overall did best at finding shelter in the first trial.

data_plot <- 
  survfit(Surv(lat_shelter, shelter_10) ~ trial,
          data = data) %>%  #filter out the accidental 6th trials
  fortify(surv.connect = TRUE) #make into a dataframe

fig_shelter_lat <- 
  
  ggplot(data_plot, aes(time, 1-surv, color = strata)) +
  
  #survival plot
  geom_step(linewidth = 1) +
  
  # axes and legends
  labs(x = "Time since start of trial (s)", y = "Proportion Frogs in Shelter") +
  
  scale_color_manual(values = rev(hcl.colors(7, palette = "BuPu")[1:5]),
                     name = "Trial number") +
  
  scale_y_continuous(expand = c(0,0), limit = c(0,0.61),
                     labels = scales::percent) +
  
  # adjust element themes
  theme_classic(base_size = 15)

fig_shelter_lat