#### Exp 2: Analysis of behavioural data #### # This experiment had 6 birds sharing two SFS units in a single room # Written by Melissa Bateson #remove any existing data/values in the global environment rm(list=ls()) # load necessary packages library(dplyr) library(ggplot2) library(reshape2) library(lubridate) library(psych) library(lme4) library(lmerTest) library(tidyr) library(PerformanceAnalytics) library(metafor) am = read.csv("Exp 2 behaviour data.csv") am$date = as.Date(am$date) # It is only valid to compare the am data between treatments because the pm data could be affected by the treatment in place: # birds might forage more actively during periods of FI but still save energy across the day. # therefore only analyse am data describe(am) # roost is the most common and most variable behaviour # roost is also by definition always visible #make a subset of just the behavioural variables b = subset(am,select=c("fly", "bathe", "perch","slide","roost","sfs","walk")) cor(b) chart.Correlation(b) # Do a PCA on the data pca.b = prcomp(b, scale = TRUE) names(pca.b) pca.b$rotation biplot(pca.b, scale = 0) std_dev = pca.b$sdev pr_var = std_dev^2 pr_var[1:7] #proportion of variance explained prop_varex = pr_var/sum(pr_var) prop_varex[1:7] # In this case the PC1 explaines 41% of the variance #scree plot plot(prop_varex, xlab = "Principal Component", ylab = "Proportion of Variance Explained", type = "b") #cumulative scree plot plot(cumsum(prop_varex), xlab = "Principal Component", ylab = "Cumulative Proportion of Variance Explained", type = "b") dim(pca.b$x) pca.out = as.data.frame(pca.b$x) # this is the PCA data am$PC1 = pca.out$PC1 # Is PC1 explained by treatment yesterday? m = lm(PC1 ~ treatment.yesterday, data = am) summary(m) # is being in shot explained by treatment? m = lm(in.shot ~ treatment.yesterday, data = am) summary(m) # birds are more likely to be in shot in FI am$active.prop = am$active/180 hist(am$active.prop) # is active explained by treatment? am$treatment.yesterday = relevel(am$treatment.yesterday, ref = "FS") m = lm(active.prop ~ treatment.yesterday, data = am) summary(m) # birds are less likely to be active in FI and FS2 #### Roosting analysis #### am$roost.prop = am$roost/180 hist(am$roost.prop) am$trans.roost = transf.arcsin(am$roost.prop) # implement the arcsin square root transformation on proportions hist(am$trans.roost) # distribution is much more normal now am$day = as.factor(am$date) # create a random effect to group data from the same day # plot a graph of proportion roosting # Plot a graph of roosting by treatment (using all days) fig = ggplot(data=am, aes(x=treatment.yesterday, y=roost.prop)) + geom_boxplot(aes(fill=treatment.yesterday)) + scale_fill_manual(values=c("cyan", "tomato", "tomato", "cyan")) + ylab("") + xlab("") + guides (fill = FALSE) + coord_cartesian(ylim = c(0,1.19)) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + annotate("segment", x=1, xend=2, y=1.06, yend=1.06, size=.5) + annotate("text", x=1.5, y=1.08, label = "NS", size=4) fig win.metafile("Exp 2 roosting.wmf") fig dev.off() # model of data am$treatment.yesterday = relevel(am$treatment.yesterday, ref = "FS") m = lmer(trans.roost ~ treatment.yesterday + (1|day), data = am) summary(m) #roosting in the morning is higher under FI; this is a behaviour that can always be seen anova(m) confint(m) mean.roost = summarise(group_by(am, treatment.yesterday), mean.roost = mean(roost.prop)) output.roost = data.frame("exp" = "Exp 2","aviary" = "single", "variable" = "prop.roost", "mean.FS" = mean.roost$mean.roost[1], "mean.FI" = mean.roost$mean.roost[2], "beta" = coef(summary(m))["treatment.yesterdayFI", "Estimate"], "se.beta" = coef(summary(m))["treatment.yesterdayFI", "Std. Error"]) # output the meta-data for roosting write.csv(output.roost,"Exp 2 roosting data for meta-analysis.csv")