--- title: "Effects previous captivity on behaviours great tits" author: "Edward Kluen" date: "07/04/2022" output: word_document editor_options: markdown: wrap: 72 --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ```{r echo=FALSE, results=FALSE, include=FALSE} setwd("") library(lme4) library (ggplot2) library (ggpubr) library (MCMCglmm) library (nadiv) library (rptR) library (tidyverse) library (broom) library(knitr) library (lattice) library(factoextra) library (corrplot) library (ggbeeswarm) library (viridis) ``` ## Short introduction This code belongs to manuscript entitled: 'Prior experience of captivity affects behavioural responses to 'novel' environments' Although we often test individuals in novel contexts to elucidate animal behaviours, these tests assume that the novelty of the test reflects underlying consistent variation among individuals in how they respond. However, information ecology theory predicts that information gathered prior to the experience from a range of sources may influence current behaviour, even if this is a different context or experience from before. Therefore, here we make use of a population of great tits (GT) of which approximately half of them were used in unrelated learning and foraging experiments. The aim of this research is to test if prior experience of handling and captivity influences three behavioural traits; exploration, breathing rate and a social mirror test. 53 great tits were captured near the Konnevesi research facilities (Central Finland) during March 2018, of which approximately half had experienced a known amount of days in captivity prior to our tests. ## Behavioural measures - Exploration - Each of the birds was subject to an exploration trial in a novel exploration room and in an exploration cage on the same day. In both arenas the number of movements in 2 minutes of the birds after they start exploring is used as the exploratory measure. - Social (mirror) interaction - For the social interaction with the mirror, we calculated the difference in time spent on each of the five perches compared to the time spent on these perches during the exploration test. We used the following metric to calculate the difference: log ((time spent on a stick during mirror exposure +1)/ (time spent on a stick during exploration +1)) where zero indicates no change and negative values indicate that the bird spent more time on the stick during the exploration period than the same mirror period. In a similar way we calculated the difference in movements during mirror versus during the exploration, we also recorded the number of pecks against the mirror and the proportion of time a bird spent looking at the mirror during the mirror test (following Carvalho et al. 2013). The proportion of time spent on each perch during the exploration trial was relative to the total amount of time a bird spent in the exploration area of the cage. For birds that were coaxed to the exploration arena of the cage, we used the five minutes after coaxing to calculate this measure. The same formula above was used to calculate difference in number of movements but instead of proportion of time on each perch, we used movements per second. A score of zero for both relative time spent on a perch and movements indicated no change in behaviour and positive values indicated that the bird spent more time on the perch, or was more active, during the mirror exposure than during exploration. - Breathing rate - We measured the time it took a bird to take 30 breaths (falls and rises of the chest) twice, using a stopwatch. Breath rate was calculated as the average of these two measurements, expressed as the number of breaths per second. The breathing rate reflects how much a bird was stressed by handling. ## PCA To reduce the multiple social response behaviour variables, we used a Principal Component Analysis (PCA). The first Principal Component (PC1; 41%) represented interaction with the mirror: number of pecks on the mirror, relative time spent on the perch closest to the mirror, and increased movement were positively loaded while relative time spent on the perch furthest from the mirror was negatively loaded (see table S3 for PC1 loadings). Subsequent PC's (PC2: 24% & PC3: 23%) were dominated by time spent looking at the mirror (PC2) and time spent on P2 and moving (PC3), however these were not used in further analyses as PC1 already included these behaviours. ```{r} #PCA exploration<- read.csv(file="datafile_explorationGT2022.csv", sep=",",dec = ".") exploration<- exploration[-c(28,32,37,84,88,93), ] # removing data of 3 birds with too short cage exploration time bird_ID: G349, G353, G358 names (exploration) #subset the data to get rid of the NA's in the bottom half of the data set. exploration_sub<- subset(exploration, test =='CAGE') # calculation of the differences in times spent looking, perching on P2 & P5 and moving between mirror and cage following the formula: # log (proportion of time spent on a perch during mirror exposure trial + C / proportion of time spent on the same perch during exploration trial + C) # where C = 1 exploration_sub$DifP2<- log((exploration_sub$P2_mir +1)/(exploration_sub$P2_expl +1)) exploration_sub$DifP5<- log((exploration_sub$P5_mir +1)/(exploration_sub$P5_expl +1)) exploration_sub$Difmove<- log((exploration_sub$mov_mir_s +1)/(exploration_sub$mov_expl_s +1)) #Making the data set for PCA including the variables: ID, looks/sec, Dif Perch2, Dif Perch5, Dif Movements and number of pecks against mirror GTmirror<- exploration_sub[, c('bird_id',"Mir_look_s","DifP2","DifP5", "Difmove", "Mir_pecks" )] # quick check summary(GTmirror) mirror.pca<-prcomp (GTmirror[,-1], scale = TRUE) summary (mirror.pca) #gbiplot(mirror.pca) # doesn't work in the current version of R....only older versions # Visualization of the variance explained by each PC fviz_eig(mirror.pca) # Results for Variables res.var <- get_pca_var(mirror.pca) # Coordinates res.var$coord # Contributions to the PCs res.var$contrib #loadings of the variables on the PC's mirror.pca$rotation ### saving the PC1 values: newdat<- data.frame(GTmirror[,1], mirror.pca$x[,1]) #the latter bit extracts the PCvalues per observation names(newdat)[1] <- "bird_id" head(newdat) length (newdat$bird_id) #write.csv(newdat, file= "PC1_values.csv") ``` ## Correlations at the phenotypic level ```{r} ##################correlations######################### #Subsetting the data by test type for exploration expl_cage<- subset(exploration, test=="CAGE") expl_room<- subset(exploration, test=="ROOM") # Preparing a dataset for the correlation test, below creating the file such that exploration in cage and room are in 2 columns (is now in 1) str(expl_room) exroom<- expl_room[,c("bird_id", "exploration_2m")] names(exroom) exroom$explor_room<- exroom$exploration_2m excage<- expl_cage[, c("bird_id", "exploration_2m", "PC1", "breath_s")] excage$explor_cage<- excage$exploration_2m exroom<- exroom[,-2] names (excage) excage<- excage[,-2] # merging the above 2 sets into 1 expl_all<- merge(exroom,excage,by="bird_id") expl_all #Phenotypic correlations cor.test(expl_all$explor_cage, expl_all$explor_room, use= "pairwise.complete.obs", method="pearson") M<-cor(expl_all[,-c(1)]) head(round(M,2),4) corrplot(M, method="number") ##getting p-values and the upper and lower Confidence intervals M1<- cor.mtest(expl_all[,-c(1)]) str (M1) head(round(M1$p,3), 4) head(round(M1$lowCI,3), 4) head(round(M1$uppCI,3), 4) ``` We test the correlations between all the behaviours at the phenotypic level and we here split the exploratory behaviour into cage and room exploration. Exploration in the room and exploration in the cage were not phenotypically correlated (Pearson's r = 0.04, CI 95%: -0.23 -- 0.31, p = 0.75, N=53), it seems that different aspects of Great tit exploration were measured in the cage and room. We found a positive (nearly significant) correlation between exploration in room and PC1 (Pearson's r = 0.26, CI 95%: -0.01 -- 0.50, p = 0.06, N = 53). Meaning that birds that explored more in the room were interacting more with the mirror. Exploration in the cage was similarly positively correlated with PC1, but weaker, where more exploratory birds in the cage were interacting more with the mirror (Pearson's r = 0.19, CI 95%: -0.08 -- 0.44, p = 0.17, N = 53). Breath rate was found not to be correlated with any of the other behaviours. ## Multi variate analyses on the 3 behaviours (results presented in Table 2 in the MS) - We test the association of the two captivity variables (and others) on each of the behavioural variables - we estimate the adjusted repeatability of the exploratory behaviour using the random effect estimates from the multivariate mixed model - We used a parameter-expanded prior that is uninformative for our model and as we have repeated measures only for the exploration variable, we set the within individual variance and within individual correlations involving behaviours with single measures to zero. We ran our analyses twice to investigate whether the length of captivity experience (i.e. continuous predictor, 'days in captivity') or the act of being held in captivity (i.e. discrete predictor, 'experience with captivity') influenced each behaviour. These models estimated the coefficients and their 95 % credible intervals of included covariates in order to assess the sign and strength of the relationship of the variables with the focal behaviour. We considered a significant effect for those covariates which credible intervals did not overlap zero. Besides captivity experience, the fixed effects structure was similar for both models and included the variables 'body mass' (scaled to their mean and unit standard deviation), 'sex' (male/female), 'age' and 'catching date' (scaled to the mean and unit standard deviation) as covariates. We also included an interaction between age and captivity history to account for variation in life experience and the possible impacts of developmental processes underlying the expressed behaviour. When estimating effects on exploration behaviour we included two additional covariates: 'test type' (cage or room) and 'test order' (whether conducted first or second). ```{r echo =T, eval=F, message = FALSE, warning = FALSE} # need to do poisson distribution for exploration because exploration is not gaussian. Apparently poisson will not be over dispersed when doing MCMCglmm # need to set uninformative priors: These are based on examples in Houslay and Wilson 2017. prior_E_Br_PC1 = list(R = list(V = diag(c(1, 0.0001, 0.0001),3,3), nu = 1.002, fix = 2), G = list(G1 = list(V = matrix(c(1,0,0, 0,1,0, 0,0,1, 0,0,0),3,3, byrow = TRUE), nu = 3, alpha.mu = rep(0,3), alpha.V = diag(c(1000, 25^2, 25^2))))) # 1000 for poisson 25^2 for gaussian # Multivarate model # where the fixed effects are: testing order and test (for the exploration trait) and bodyweight, days in captivity and sex for all behavioural traits # we use a poisson distribution for exploration and gaussian for all others. 3.050.000 iterations, 50.000 burn in and 3000 thinning # This gives that the coef. estimates are based on 1000 calculations. #(Running the following code may last over half an hour) mcmc_E_Br_PC1_capt <- MCMCglmm(cbind(exploration_2m, scale (breath_s), PC1) ~ trait-1 + at.level(trait,1):order_c.r+ at.level(trait,1):test + trait:scale(bodyweight) + trait:scale(days_captivity)+ trait: scale(March_day) +trait:sex + trait:age + trait:scale(days_captivity):age, random =~ us(trait):bird_id, rcov =~ us (trait):units, family = c("poisson","gaussian", "gaussian"), prior = prior_E_Br_PC1, nitt=3050000, burnin=50000, thin=3000, verbose = TRUE, pr = TRUE, #This saves the posterior distribution of the individual random effects takes up space data = as.data.frame(exploration)) # to see the posterior distributions of the model: plot(mcmc_E_Br_PC1_capt) plot (mcmc_E_Br_PC1_capt$VCV) #trace and density plots of the (co)variances of the random effects # to get at the coefficients of the random and fixed effects summary(mcmc_E_Br_PC1_capt) # Similar model, now with captivity experience expressed as y/n mcmc_E_Br_PC1_exp <- MCMCglmm(cbind(exploration_2m, breath_s, PC1) ~ trait-1 + at.level(trait,1):test + at.level(trait,1):order_c.r + trait:scale(bodyweight) + trait:experience+ trait:sex + trait:scale(March_day)+trait:age + trait:experience:age, random =~ us(trait):bird_id, rcov =~ us (trait):units, family = c("poisson","gaussian", "gaussian"), prior = prior_E_Br_PC1, #same prior settings as above MMM nitt=3050000, burnin=50000, thin=3000, verbose = TRUE, pr = TRUE, #This saves the posterior distribution of the individual random effects takes up space data = as.data.frame(exploration)) # to see the posterior distributions of the model: plot (mcmc_E_Br_PC1_exp) #to check if ran nicely plot(mcmc_E_Br_PC1_exp$VCV) #trace and density plots of the (co)variances of the random effects # to get at the coefficients of the random and fixed effects summary(mcmc_E_Br_PC1_exp ) ``` To test whether the association of captivity experience with exploratory behaviour is dependent on the arena-type (cage/room) in which exploration was tested. We checked whether test type influenced effects of captivity history (assessed using an interaction term in a uni-variate model; results in table S4 in Supplement) but as it was not significant, we did not include this interaction in the MMMs. ```{r} # 2 Univariate mixed models to test the interaction term of captivity experience (both continuous and discrete) and arena type ('test' (room/cage)) # These results are presented in Table S4 in the supplement of the MS. prior_Ex <- list(R = list(V = 1, nu = 0.002), G = list(G1 = list(V =1, nu = 0.002, alpha.mu = 0, alpha.V = 1000))) #continuous mcmc_E_capt <- MCMCglmm(exploration_2m ~ order_c.r+ test + scale(days_captivity):test + scale(bodyweight) + scale(days_captivity)+ scale (March_day) + sex + age + scale(days_captivity):age, random =~ bird_id, family = "poisson", prior = prior_Ex, nitt=3050000, burnin=50000, thin=3000, verbose = TRUE, pr = TRUE, data = as.data.frame(exploration)) plot(mcmc_E_capt) summary(mcmc_E_capt) #discrete mcmc_E_experience <- MCMCglmm(exploration_2m ~ order_c.r+ test + experience:test + scale(bodyweight) + experience+ scale (March_day) +sex + age + experience:age, random =~ bird_id, family = "poisson", prior = prior_Ex, nitt=3050000, burnin=50000, thin=3000, verbose = TRUE, pr = TRUE, data = as.data.frame(exploration)) plot(mcmc_E_experience) summary(mcmc_E_experience) #in both models the interaction between captivity history and arena type is not significant. ``` ## Does experience with captivity explain any of the variation in the behavioural variables? - Experience with captivity expressed as *days in captivity* or by the *experience with captivity (y/n)* is positively associated with the exploratory behaviour. Meaning that birds that have been longer in captivity or had captive experience were also more explorative. - the interaction of *days in captivity* and *age* is significantly associated with exploratory behaviour. Where older birds that with more captive days are less exploratory than younger birds with similar captive experience, this interaction is also significant for the interaction of *experience with captivity* and *age*. - we do not find that prior experience of captivity is associated with the *social interaction with the mirror (PC1)* or breathing rate (handling stress). ## Some pre-work for the graphs ```{r} exploration$age<- as.factor(exploration$age) levels(exploration$age) <- c("young","old") exploration$experience<- as.factor (exploration$experience) levels(exploration$experience) <- c("no","yes") # To be able to plot the mean and SE's for each of the plots I calculate these values (tidyverse) to plot into the violin plots # because we used a poisson distribution in the MMModels I will plot the exploration here log transformed explSummary <- exploration %>% group_by(experience) %>% summarize(expl_mean = mean(log(exploration_2m+1)), expl_se = sqrt(var(log(exploration_2m+1))/length(exploration_2m))) explSummary explSum_age <- exploration %>% group_by(experience, age) %>% summarize(expl_mean = mean(log(exploration_2m+1)), expl_se = sqrt(var(log(exploration_2m+1))/length(exploration_2m+1))) explSum_age #Subset the data to get rid of the NA's in the dataset and then calculate the SE's # each of the following summaries contains the mean value of a behaviour by experience and age. expl_cage<- subset(exploration, test == "CAGE") #breathing rate brSummary <- expl_cage %>% group_by(experience) %>% summarize(br_mean = mean(breath_s), br_se = sqrt(var(breath_s)/length(breath_s))) brSummary #PC1 pc1Summary <- expl_cage %>% group_by(experience) %>% summarize(pc1_mean = mean(PC1), pc1_se = sqrt(var(PC1)/length(PC1))) pc1Summary ``` ## Graphs (Figure 2 in the MS) of the relationship of each behaviour and captivity (length) and experience. ```{r fig.cap="The relationship of each behaviour with days in captivity and experience with captivity where data are plotted by age ", echo =FALSE, warning = FALSE, message = FALSE, eval=TRUE, fig.height = 18, fig.width = 12} #Plots #A) #sign. Interaction of days in captivity and age on Exploration (lines), with overall days captivity-exploration relationship (black line) capt3<- ggplot(exploration, aes (x= scale(days_captivity), y= log(explor_2m+1))) + geom_point(aes(shape = test, color = age), size = 2.5, position = "jitter") + scale_color_manual(values=c("#F5793A","#85C0F9"))+ geom_smooth(aes(group = age, color = age), method = 'lm', se = TRUE, alpha = 0.15, lty= "dashed")+ geom_smooth(aes(x= scale(days_captivity), y= log(explor_2m+1)), method = "lm", se=FALSE, linetype="solid", color="darkgrey")+ theme_classic(base_size = 10)+ xlim(-1, 3)+ xlab("Days in captivity (scaled)") + ylab("log Exploration (moves 2min)") #B) # relationship of experience with captivity on Exploration by age capt4<- ggplot(exploration, aes (x= experience, y= log(explor_2m+1), fill = age)) + geom_violin(trim=TRUE, alpha = 0.4, position="dodge", show.legend= FALSE)+ scale_fill_manual(values=c("light grey","dark grey"))+ scale_x_discrete(limits=c("no","yes"))+ geom_quasirandom(aes(x=experience, y= log(explor_2m+1), group= age, color =age, shape = test) , size =2.5, dodge.width = 0.9, width = 0.15, varwidth = TRUE)+ scale_color_manual(values=c("#F5793A","#85C0F9"))+ geom_point(aes(y = expl_mean, group = age), color = "black", cex= 1.5, position = position_dodge(width = 0.9), data = explSum_age) + geom_errorbar(aes(y = expl_mean, ymin = expl_mean-expl_se, ymax = expl_mean+expl_se, group = age), color = "black", cex = 0.2, width = 0.2, position = position_dodge(width = 0.9), data = explSum_age)+ theme_classic(base_size = 10) + xlab("Experience of captivity") + ylab("log Exploration (moves 2min)") #C) relationship of days in captivity with PC1, by age capt9<- ggplot(exploration, aes (x= scale(days_captivity), y= PC1)) + geom_point(shape = 18, color = "darkgrey",size = 2.5, position = "jitter") + scale_color_manual(values=c("#F5793A","#85C0F9"))+ geom_smooth(method = 'lm', se = FALSE, color= "darkgrey", lty= 9)+ theme_classic(base_size = 12)+ xlim(-1, 3)+ xlab("days in captivity (scaled)") + ylab("PC1 (mirror interaction)") #D) relationship of experience and PC1, not sign. interaction with age capt12<- ggplot(exploration, aes (x= experience, y= PC1, fill = experience)) + geom_violin(trim=TRUE, alpha = 0.4, position="dodge", show.legend= FALSE)+ scale_x_discrete(limits=c("no","yes")) + scale_fill_manual(values=c("light grey","dark grey"))+ geom_quasirandom(aes(x=experience, y= PC1), size =2.5, shape=18, color = "darkgrey", dodge.width = 1, width = 0.2, varwidth = TRUE)+ geom_point(aes(y = pc1_mean), color = "black", size = 1.8, position = position_dodge(width = 0.9),show.legend= FALSE, data = pc1Summary) + geom_errorbar(aes(y = pc1_mean, ymin = pc1_mean-pc1_se, ymax = pc1_mean+pc1_se), color = "black", width = 0.2,position = position_dodge(width = 0.9), show.legend= FALSE,data = pc1Summary)+ theme_classic(base_size = 12)+ xlab("Experience of captivity") + ylab("PC1 (mirror interaction)") #E)relationship of days in captivity with breath rate, by age capt5<- ggplot(exploration, aes (x= scale(days_captivity), y= breath_s)) + geom_point(shape = 18, color = "darkgrey", size = 2.5, position = "jitter") + scale_color_manual(values=c("#F5793A","#85C0F9"))+ geom_smooth(method = 'lm', se = FALSE, color = "darkgrey", lty= 9)+ theme_classic(base_size = 12)+ xlim(-1, 3)+ xlab("days in captivity (scaled)") + ylab("Breathing rate (breath/s)") # F)relationship of experience with captivity with breath rate by age capt8<- ggplot(exploration, aes (x= experience, y= breath_s, fill = experience)) + geom_violin(trim=TRUE, alpha = 0.4, position = "dodge", color = "black", show.legend= FALSE)+ scale_fill_manual(values=c("light grey","dark grey"))+ scale_x_discrete(limits=c("no","yes")) + geom_quasirandom(aes(x=experience, y= breath_s), size =2.5, shape=18, color = "darkgrey", dodge.width = 1, width = 0.2, varwidth = TRUE)+ geom_point(aes(y = br_mean), color = "black", size = 1.8, position = position_dodge(width = 0.9),show.legend= FALSE, data = brSummary) + geom_errorbar(aes(y = br_mean, ymin = br_mean-br_se, ymax = br_mean+br_se), color = "black", width = 0.2,position = position_dodge(width = 0.9),show.legend= FALSE, data = brSummary)+ theme_classic(base_size = 12)+ xlab("Experience of captivity") + ylab("Breathing rate (breath/s)") #Making the plot all in one ggarrange(capt3, capt4, capt9, capt12, capt5, capt8, labels = c("A", "B", "C","D" ,"E" ,"F"),common.legend = TRUE, ncol = 2, nrow = 3,legend = "bottom", label.x = 0.0, label.y = 1) ``` ## Repeatability of exploration (Table 2) ```{r echo =T, eval=F, message = FALSE, warning = FALSE} #Repeatability of Exploratory behaviour mcmc_prop_E<- mcmc_E_Br_PC1_capt $VCV[,"traitexploration_2m:traitexploration_2m.bird_id"]/ (mcmc_E_Br_PC1_capt$VCV[,"traitexploration_2m:traitexploration_2m.bird_id"] + mcmc_E_Br_PC1_capt$VCV[,"traitexploration_2m:traitexploration_2m.units"]) plot(mcmc_prop_E) #repeatability mean(mcmc_prop_E) #Credible intervals HPDinterval(mcmc_prop_E) ``` # Reduction of repeatability due the inclusion of captivity experience (Figure S3 in the supplement of the MS) Below we show how much the repeatability of exploration changes when we take into account captivity history. Of course we have just seen that repeatability is not sigificant (coming from the MMM), we use here the package rptR to test the adj repeatability of exploration with and without the captive history (days_captivity). ```{r echo =T, eval=F, message = FALSE, warning = FALSE} #Repeatability exploration behaviour, using package rptR names (exploration) rep_expl <- rpt(exploration_2m~ test+ order_c.r+ (1 | bird_id), grname = c("bird_id", "Residual"), data = exploration, datatype = "Poisson", nboot = 1000, npermut = 0, ratio = TRUE) print(rep_expl) #Repeatability exploration with day_captivity in the model rep_expl_d_act <- rpt(exploration_2m~ test+ order_c.r+ days_captivity+ (1 | bird_id), grname = c("bird_id", "Residual"), data = exploration, datatype = "Poisson", nboot = 1000, npermut = 0, ratio = TRUE) print(rep_expl_d_act) par(mfrow= c(1,2)) p1<- plot(rep_expl, cex.main = 1) mysubtitle = expression(italic("adjusted repeatability without captive experience, R = 0.22, CI = [0, 0.53], P = 0.085 [LRT] ")) mtext(side=3, line=0, cex=0.8, mysubtitle) p2<- plot(rep_expl_d_act, cex.main = 1) mysubtitle1 = expression(italic("adjusted repeatability with captive experience, R = 0.198, CI = [0, 0.48], P = 0.11 [LRT] ")) mtext(side=3, line=0, cex=0.8, mysubtitle1) ```