#### Expt 3: Analysis of starling guano samples ####
# Written by Melissa Bateson
#remove any existing data/values in the global environment
rm(list=ls())
# load necessary packages
library(dplyr)
library(ggplot2)
library(lme4)
library(lmerTest)
library(reshape2)
# read in and merge the data
d = read.csv("Exp 3 bomb calorimeter data.csv")
#### Reliability analysis of technical replicates to examine measurement error ####
new = data.frame(d$date, d$tech.rep, d$energy)
new = dcast(data = new, d.date~ ...) # make the dataframe wide
plot(new$`1`, new$`2`)
cor.test(new$`1`, new$`2`)
# correlation between technical replicates is 0.046, but little variation therefore hard to interpret
#### Models and graphs of effects of FI treatment ####
d$treatment = relevel(d$treatment, ref = "FS1") # change ref level to get phases in the right order
#d$treatment.wk = relevel(d$treatment.wk, ref = "FS1")
# calculate the percentage decline in energy density under FI using the means of the raw data
summarise(group_by(d,treatment), mean.energy.density = mean(energy), sd.energy.density = sd(energy))
percent.decline = ((15.3-15.2)/15.3)*100 # 0.65% decline in the energy density of faeces under FI
# calculate cohen's d for effect of FI
cohensd = (15.3-15.2)/ sqrt(((0.223^2 + 0.181^2)/2)) # 0.49 i.e. a medium effect size
# Is there a significant effect of phase on energy?
# need to add a random effect of sample to account for technical replicates
d$sample = as.factor(d$date)
d$treatment = relevel(d$treatment, ref = "FS1")
m1 = lmer(energy ~ treatment + (1|sample), data=d)
summary(m1) # sig diff between FS1 and FS2
confint(m1)
anova(m1) # for overall effect of phase
d$treatment = relevel(d$treatment, ref = "FI") # change ref level to get contrast between phases 2 and 3
m1b = lmer(energy ~ treatment + (1|sample), data=d)
summary(m1b) # sig diff between FS1 and FI and FS2
confint(m1b)
# Plot energy by each main phase of the expt
d$treatment.wk = relevel(d$treatment.wk, ref = "FS1") # change ref level to get phases in the right order
fig1 = ggplot(data=d, aes(x=as.factor(treatment.wk), y=energy)) +
geom_boxplot(aes(fill=treatment.wk)) +
scale_fill_manual(values=c("cyan", "tomato", "cyan", "cyan")) +
guides(fill=FALSE) +
ylim(14.5,15.85) +
ylab("Energy density of guano (MJ/kg)") +
xlab("") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
annotate("segment", x=2, xend=3.5, y=15.75, yend=15.75, size=1) +
annotate("text", x=2.75, y=15.77, label = "*", size=4) +
annotate("segment", x=1, xend=3.5, y=15.80, yend=15.80, size=1) +
annotate("text", x=2.25, y=15.82, label = "**", size=5)
fig1
win.metafile("Exp 3 guano energy by treatment.wmf")
fig1
dev.off()
# Is there a decline in energy as FI progresses?
d1 = d[d$treatment!="FS1",] # take subset of days with FI and FS2
d1$day.FI[d1$treatment=="FI"] = d1$day.treat[d1$treatment=="FI"]
d1$day.FI[d1$treatment!="FI"] = 7+d1$day.treat[d1$treatment!="FI"]
m2 = lmer(energy ~ day.FI + (1|sample), data=d1)
summary(m2)
anova(m2) # Significantly decline in energy with increasing days following the start of FI
confint(m2) # calculate the 95% CI for the slope
fig2 = ggplot(data=d1, aes(x=day.FI, y=energy)) +
geom_point(aes(colour = treatment), size = 2.5) +
geom_smooth(method = "lm", colour = "black", alpha = 1) +
labs(colour = "Treatment") +
ylim(14.5,15.75) +
ylab("Energy density of guano (MJ/kg)") +
xlab("Day since start of FI treatment") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = c(.08,.13)) +
theme(legend.background = element_rect(colour = "black"))
fig2
win.metafile("Exp 3 guano energy by day FI.wmf")
fig2
dev.off()
# Calculate the energy density of food
food = f1[f1$treatment =="food",]
food.energy = mean(f1$energy) # mean energy content of food.