---
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)
```