#### SFS Expt 1: 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 the data
d = read.csv("Exp 1 bomb calorimeter data.csv")
#### Reliability analysis of technical replicates to examine measurement error ####
new = data.frame(d$sample, d$tech.rep, d$energy)
new = dcast(data = new, d.sample~ ...) # make the dataframe wide
plot(new$`1`, new$`2`)
cor.test(new$`1`, new$`2`)
# correlation between technical replicates is only 0.37, 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.1)/15.3)*100 # 1.31% decline in the energy density of faeces under FI
# calculate cohen's d for effect of FI
cohensd = (15.3-15.1)/ sqrt(((0.249^2 + 0.195^2)/2)) # 0.89 i.e. a large effect size
# Is there a significant effect of phase on energy?
# need to include random effects for aviary and technical replicates nested within aviary
m1 = lmer(energy ~ treatment + (1|room/sample), data=d)
summary(m1) # sig diff between FS1 and FI 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|room/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 = relevel(d$treatment, ref = "FS1") # change ref level to get phases in the right order
sum.tech2 = summarise(group_by(d,treatment, treatment.wk,room), mean.energy = mean(energy))
# now calculate overall mean for each room and merge this back into the above dataframe
mean.room = summarise(group_by(sum.tech2,room), mean.room.energy = mean(mean.energy))
sum.tech2 = merge(sum.tech2,mean.room, by="room")
# within-room centering followed by adding back the grant mean so that phase effects are relocated relative to the grand mean
sum.tech2$norm.energy = (sum.tech2$mean.energy - sum.tech2$mean.room.energy) + mean(sum.tech2$mean.energy)
fig1 = ggplot(data=sum.tech2, aes(x=as.factor(treatment.wk), y=norm.energy)) +
geom_boxplot(aes(fill=treatment.wk)) +
scale_fill_manual(values=c("cyan", "tomato", "tomato", "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=1, xend=2.5, y=15.6, yend=15.6, size=1) +
annotate("text", x=1.75, y=15.62, label = "*", size=5) +
annotate("segment", x=1, xend=4, y=15.67, yend=15.67, size=1) +
annotate("text", x=2.5, y=15.69, label = "**", size=5)
fig1
win.metafile("Exp 1 guano energy by treatment.wmf")
fig1
dev.off()
# Is there a decline in energy as FI progresses?
d1 = d[d$treatment=="FI",] # take subset of days with FI
m2 = lmer(energy ~ day.treat + (1|room/sample), data=d1)
summary(m2)
anova(m2) # Significantly decline in energy with increasing days of FI
confint(m2) # calculate the 95% CI for the slope
m3 = lm(energy ~ day.treat, data=d1)
summary(m3)
d2 = d[d$treatment!="FS1",] # take subset of days with FI
d2$day.FI = d2$day.expt-7
m2 = lmer(energy ~ day.FI + (1|room/sample), data=d2)
summary(m2)
anova(m2) # Marginally non-significantly decline in energy with increasing days of FI
fig2 = ggplot(data=d1, aes(x=day.treat, y=energy)) +
geom_point(colour = "coral", size = 2.5) +
geom_smooth(method = "lm", colour = "black") +
guides(shape=FALSE) +
labs(shape = "Room") +
ylab("Energy density of faeces (MJ/kg)") +
xlab("Day of FI treatment") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig2
fig3 = ggplot(data=d2, 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("Days since start of FI treatment") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = c(.08,.12)) +
theme(legend.background = element_rect(colour = "black"))
fig3
win.metafile("Exp 1 guano energy by day FI.wmf")
fig3
dev.off()