#### 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")