#### Exp 4: Analysis of food eaten and masses #### # This experiment had 6 birds housed in 3 pairs each sharing a single SFS unit # Written by Melissa Bateson #remove any existing data/values in the global environment rm(list=ls()) library(dplyr) library(ggplot2) library(reshape2) library(lubridate) library(psych) library(lme4) library(lmerTest) library(tidyr) # Useful functions # create a new version of length which can handle NA's: if na.rm==T, don't count them length2 <- function (x, na.rm=FALSE) { if (na.rm) sum(!is.na(x)) else length(x) } # function to turn hr:min:sec to minutes since midnight as a decimal decimateTime=function(time) { time=as.numeric(unlist(strsplit(time, ":"))) time = time[1]*60+time[2]+time[3]/60 return(time) } #### Anaysis of daily food consumption #### f1 = read.csv("Exp 4 food consumption data.csv") f1$date = as.Date(f1$date) f1$total.eaten = f1$total_eaten f1$total_eaten = NULL start = c("2018-03-12","2018-04-25","2018-05-15") # the dates of catching for each pair start1 = c("2018-03-26","2018-05-12","2018-06-10") # start of p(r) = 1 (changeover date) start0.4 = c("2018-04-07","2018-05-21","2018-06-21") # start of p(r) = 0.4 (changeover date) start0.2 = c("2018-04-17","2018-05-29","2018-06-30") # start of p(r) = 0.2 (changeover date) pair = c(1,2,3) start.date = data.frame(pair,start,start1,start0.4,start0.2) start.date$start = as.Date(start.date$start) start.date$start1 = as.Date(start.date$start1) start.date$start0.4 = as.Date(start.date$start0.4) start.date$start0.2 = as.Date(start.date$start0.2) f1 = merge(f1,start.date, by = "pair") f1$day.exp = f1$date-f1$start # the number of days in the expt at the time of a mass # the number of days in the current treatment at the time of a mass f1$day.treat[f1$probability=="1"] = (f1$date-f1$start1)[f1$probability=="1"] f1$day.treat[f1$probability=="0.4"] = (f1$date-f1$start0.4)[f1$probability=="0.4"] f1$day.treat[f1$probability=="0.2"] = (f1$date-f1$start0.2)[f1$probability=="0.2"] # code the treatments f1$treatment[f1$probability==1] = "FS" f1$treatment[f1$probability==0.4] = "FIlow" f1$treatment[f1$probability==0.2] = "FIhigh" f1$treatment = factor(f1$treatment, levels = c("FS", "FIlow", "FIhigh")) # Basic stats per bird mean.eaten.pair = summarise(group_by(f1, treatment, pair), mean.food.eaten = mean(total.eaten/2, na.rm = TRUE)) mean.eaten = summarise(group_by(mean.eaten.pair, treatment), mean.eaten.bird = mean(mean.food.eaten), sd.food.eaten = sd(mean.food.eaten), n.food.eaten = length(mean.food.eaten)) # calculate the percentage drop in consumption between FI and FS1 percent.drop.FIlow = 100*(mean.eaten$mean.eaten.bird[1]-mean.eaten$mean.eaten.bird[2])/mean.eaten$mean.eaten.bird[1] percent.drop.FIhigh = 100*(mean.eaten$mean.eaten.bird[1]-mean.eaten$mean.eaten.bird[3])/mean.eaten$mean.eaten.bird[1] # within-based pair centering based on raw values for plotting mean.pair = summarise(group_by(f1, pair), mean.food.pair = mean((total.eaten/2), na.rm=T)) f1 = merge(f1,mean.pair, by="pair") # within-pair centering followed by adding back the grand mean so that treatment effects are relocated relative to the grand mean f1$norm.eaten = ((f1$total.eaten/2) - f1$mean.food.pair)+(mean((f1$total.eaten), na.rm=T)/2) f1$norm.eaten = ((f1$total.eaten/2) - f1$mean.food.pair)+mean(mean.eaten$mean.eaten.bird) # Plot a box plot of consumption by treatment (using all days) f1$treatment = relevel(f1$treatment, ref = "FS") fig = ggplot(data=f1, aes(x=treatment, y=norm.eaten)) + geom_boxplot(aes(fill=treatment)) + scale_fill_manual(values=c("cyan", "tomato", "tomato")) + #ylab("Food consumption (g/day)") + ylab("") + xlab("") + ylim(9,22) + scale_x_discrete(labels = c("FS" = "FS", "FIlow"=expression(paste(FI[low])), "FIhigh"=expression(paste(FI[high])))) + guides (fill = FALSE) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + annotate("segment", x=1, xend=2, y=20, yend=20, size=.5) + annotate("text", x=1.5, y=20.2, label = "NS", size=3) + annotate("segment", x=2, xend=3, y=19, yend=19, size=.5) + annotate("text", x=2.5, y=19.2, label = "**", size=4) + annotate("segment", x=1, xend=3, y=21, yend=21, size=.5) + annotate("text", x=2, y=21.2, label = "***", size=4) fig win.metafile("Exp 4 consumption.wmf") fig dev.off() # Model amount eaten f1$treatment = relevel(f1$treatment, ref = "FS") m = lmer((total.eaten/2) ~ treatment + (1|pair), data = f1) summary(m) # The birds eat significantly less under FIhigh compared to FS confint(m) anova(m) m2 = lmer(scale(total.eaten/2) ~ treatment + (1|pair), data = f1) summary(m2) f1$treatment = relevel(f1$treatment, ref = "FIlow") # change ref level to get comparison of two levels of FI m1 = lmer((total.eaten/2) ~ treatment + (1|pair), data = f1) summary(m1) # The birds eat significantly less under FIhigh compared to FIlow confint(m1) # output the meta-data output.food = data.frame("exp" = "Exp 4","aviary" = "all", "variable" = "food.eaten", "mean.FS" = mean.eaten$mean.eaten.bird[1], "mean.FI" = mean.eaten$mean.eaten.bird[3], "beta" = coef(summary(m))["treatmentFIhigh", "Estimate"], "se.beta" = coef(summary(m))["treatmentFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2))["treatmentFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2))["treatmentFIhigh", "Std. Error"]) # Now calculate the effects of FI on consumption by aviary for the by aviary meta-analysis f1$treatment = relevel(f1$treatment, ref = "FS") # change ref level # pair 1 f1.1 = f1[f1$pair==1, ] m.1 = lm(total.eaten/2 ~ treatment, data=f1.1) summary(m.1) m2.1 = lm(scale(total.eaten/2) ~ treatment, data=f1.1) summary(m2.1) output.eaten.1 = data.frame("exp" = "Exp 4","aviary" = "aviary 1", "variable" = "food.eaten", "mean.FS" = mean.eaten.pair$mean.food.eaten[1], "mean.FI" = mean.eaten.pair$mean.food.eaten[7], "beta" = coef(summary(m.1))["treatmentFIhigh", "Estimate"], "se.beta" = coef(summary(m.1))["treatmentFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.1))["treatmentFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.1))["treatmentFIhigh", "Std. Error"]) # pair 2 f1.2 = f1[f1$pair==2, ] m.2 = lm(total.eaten/2 ~ treatment, data=f1.2) summary(m.2) m2.2 = lm(scale(total.eaten/2) ~ treatment, data=f1.2) summary(m2.2) output.eaten.2 = data.frame("exp" = "Exp 4","aviary" = "aviary 2", "variable" = "food.eaten", "mean.FS" = mean.eaten.pair$mean.food.eaten[2], "mean.FI" = mean.eaten.pair$mean.food.eaten[8], "beta" = coef(summary(m.2))["treatmentFIhigh", "Estimate"], "se.beta" = coef(summary(m.2))["treatmentFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.2))["treatmentFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.2))["treatmentFIhigh", "Std. Error"]) # pair 3 f1.3 = f1[f1$pair==3, ] m.3 = lm(total.eaten/2 ~ treatment, data=f1.3) summary(m.3) m2.3 = lm(scale(total.eaten/2) ~ treatment, data=f1.3) summary(m2.3) output.eaten.3 = data.frame("exp" = "Exp 4","aviary" = "aviary 3", "variable" = "food.eaten", "mean.FS" = mean.eaten.pair$mean.food.eaten[3], "mean.FI" = mean.eaten.pair$mean.food.eaten[9], "beta" = coef(summary(m.3))["treatmentFIhigh", "Estimate"], "se.beta" = coef(summary(m.3))["treatmentFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.3))["treatmentFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.3))["treatmentFIhigh", "Std. Error"]) #### Mass analysis #### #Read in .csv file d = read.csv("Exp 4 mass data.csv") d$bird = d$ring levels(d$bird) d$ring = NULL d$room = as.factor(d$pair) #Specify variable type d$balance<-factor(d$balance) d$RFID<-factor(d$RFID) d$hopper_status<-factor(d$hopper_status) d$pair<-factor(d$pair) d$pair.id<-factor(d$pair.id) # code the treatments d$treatment[d$probability==1] = "FS" d$treatment[d$probability==0.4] = "FIlow" d$treatment[d$probability==0.2] = "FIhigh" d$treatment = factor(d$treatment, levels = c("FS", "FIlow", "FIhigh")) #Create a mass variable in g d$mass_g<-d$mass_kg*1000 #Create a datetime variable d$date_time = ymd_hms(paste(d$date, d$time)) #Specify date variable d$date<-as.Date(d$date,format="%Y-%m-%d") #Create an hour variable that allows us to select masses according to the hour they were recorded d$bin.1h<-as.numeric(format(d$date_time,"%H")) # Added by MB: create a decimal time variable (minutes since midnight) d$time.char = as.character(d$time) d$time.dec = (as.numeric(as.POSIXct(paste("2014-01-01", d$time.char))) - as.numeric(as.POSIXct("2014-01-01 0:0:0")))/60 # this code converts h:m:s to decimal minutes from midnight #Create a 'day.data' variable and merge it into the main dataframe d.date<-unique(subset(d,select=c("date","pair","probability"))) d.date<-d.date%>% group_by(pair,probability)%>% mutate(day.data=seq_along(date)) start = c("2018-03-12","2018-04-25","2018-05-15") # the dates of catching for each pair start1 = c("2018-03-26","2018-05-12","2018-06-10") # start of p(r) = 1 (changeover date) start0.4 = c("2018-04-07","2018-05-21","2018-06-21") # start of p(r) = 0.4 (changeover date) start0.2 = c("2018-04-17","2018-05-29","2018-06-30") # start of p(r) = 0.2 (changeover date) pair = c(1,2,3) start.date = data.frame(pair,start,start1,start0.4,start0.2) start.date$start = as.Date(start.date$start) start.date$start1 = as.Date(start.date$start1) start.date$start0.4 = as.Date(start.date$start0.4) start.date$start0.2 = as.Date(start.date$start0.2) d = merge(d,start.date, by = "pair") d.date = subset(d.date, select = c("date","day.data")) d = merge(d,d.date, by = "date") d$day.exp = d$date-d$start # the number of days in the expt at the time of a mass # the number of days in the current treatment at the time of a mass d$day.treat[d$probability=="1"] = (d$date-d$start1)[d$probability=="1"] d$day.treat[d$probability=="0.4"] = (d$date-d$start0.4)[d$probability=="0.4"] d$day.treat[d$probability=="0.2"] = (d$date-d$start0.2)[d$probability=="0.2"] #Exclude any times between 11:00:00 - 12:00:00 as this was when husbandry occurred, so #any masses recorded here may be unreliable as birds a) had a lib mealworms and so were possibly fatter #than normal and b) birds probably didn't visit the feeder normally anyway, as were worried about #experimenter presence d = subset(d,bin.1h!=11.0) #Note, that the lights come on between 06:00-21:00 and the birds' masses can be recorded by starfeeder #throughout those times. Thus, it means that we can bin our masses by hour. However, notice #that for pecks, the continuous program starts at 06:15 when lights are on fully and ends #at 20:45 when the lights start to dim, so we can't use hour in the same way if we want an equal number of bins #The continuous program behaved in this way to not disturb the birds during dawn/dusk # number of weights before 06:15 and after 20:45 summarise(group_by((subset(d,time.dec<375)),pair.id,treatment), n = length2(time.dec, na.rm = TRUE)) summarise(group_by((subset(d,time.dec<375)),treatment), n = length2(time.dec, na.rm = TRUE)) summarise(group_by((subset(d, time.dec>1245)),pair.id,treatment), n = length2(time.dec, na.rm = TRUE)) summarise(group_by((subset(d, time.dec>1245)),treatment), n = length2(time.dec, na.rm = TRUE)) # no. masses with data excluded d = subset(d, (time.dec>375&time.dec<1245)) # excluding masses when SFS not delivering food summarise(group_by(d,pair.id,treatment), n = length2(time.dec, na.rm = TRUE)) summarise(group_by(d,treatment), n = length2(time.dec, na.rm = TRUE)) no.masses = summarise(group_by(d,bird), n = length2(mass_g, na.rm = TRUE)) mean(no.masses$n)/21 # mean number of masses/bird/day = 91.91 mean(no.masses$n)/(21*14.5) # mean number of masses/bird/hr = 6.34 # create a unique column with a unique row number (to be used later when remerging dataframes) d$row.number = seq(1,length(d$date)) #### Polynomial fitting #### # The aim of this part of the script is to fit the mass data obtained for each individual bird on each day of the experiment. # The reason for this is to obtain a model of the mass change over the day that eliminates the noise in the measurements. # We can then use the fitted values from the model in subsequent analyses (e.g. the maximum fitted mass in a day; the fitted mass at 1800 hrs) # set the mass band either side of the fitted function for retaining masses g.filter = 3 # I have finally settled on quite a high value that only excludes very extreme masses. # make a new data frame to accumulate cleaned and fitted masses in d.clean = data.frame(time.dec = numeric(0),time_s.clean = numeric(0), mass_g.clean= numeric(0), m.fit = numeric(0)) #make empty vectors for results room = NULL bird = NULL treatment = NULL day.exp = NULL day.data = NULL day.treat = NULL no.masses = NULL no.masses.clean = NULL linear.p = NULL quadratic.p = NULL cubic.p = NULL r_sq.fit = NULL AIC.fit = NULL mass.earliest = NULL mass.latest = NULL mass.max = NULL mass.min = NULL mass.dawn = NULL mass.7 = NULL mass.8 = NULL mass.9 = NULL mass.10 = NULL mass.11 = NULL mass.12 = NULL mass.13 = NULL mass.14 = NULL mass.15 = NULL mass.16 = NULL mass.17 = NULL mass.18 = NULL mass.19 = NULL mass.20 = NULL mass.dusk = NULL am.gain = NULL time.max.mass = NULL time.min.mass = NULL time.earliest.mass = NULL time.latest.mass = NULL row = 1 for(i in c("BG/RR", "BG/YY", "BY/BB", "BY/PP", "GR/BB", "GR/PP")){ print(i) for(j in c("FS","FIlow", "FIhigh")){ print(j) for(k in 1:7){ print(k) bird[row] = i treatment[row] = j day.data[row] = k s = subset(d, (bird == i & treatment == j & day.data == k)) room[row] = as.character(s$room[1]) day.exp[row] = as.numeric(s$day.exp[1]) day.treat[row] = s$day.treat[1] no.masses[row] = length(s$mass_g) # save the number of masses that day # check to see that there is enough data to do the fitting # 'degree' must be less than number of unique points; need degree+1 points for fit that is not just the data # don't fit the data if there are no masses in the first two hours of foraging (495 mins) if ((length(s$mass_g)<10)|(min(s$time.dec)>495)) { # Set fit parameters to NA linear.p[row] = NA quadratic.p[row] = NA cubic.p[row] = NA r_sq.fit[row] = NA AIC.fit[row] = NA mass.earliest[row] = NA mass.latest[row] = NA mass.max[row] = NA mass.min[row] = NA mass.dawn[row] = NA mass.7[row] = NA mass.8[row] = NA mass.9[row] = NA mass.10[row] = NA mass.11[row] = NA mass.12[row] = NA mass.13[row] = NA mass.14[row] = NA mass.15[row] = NA mass.16[row] = NA mass.17[row] = NA mass.18[row] = NA mass.19[row] = NA mass.20[row] = NA mass.dusk[row] = NA am.gain[row] = NA time.max.mass[row] = NA time.min.mass[row] = NA time.earliest.mass[row] = NA time.latest.mass[row] = NA row = row + 1 } if ((length(s$mass_g) >=10)&(min(s$time.dec)<=495)) { m = lm(mass_g ~ poly(time.dec,degree=3), data = s) # NB this fits a cubic orthogonal polynomial s$m.fit.raw = predict(m) # calculate a fixed mass band for inclusion either side of the fitted slope s$m.lb=(s$m.fit.raw - g.filter) s$m.ub=(s$m.fit.raw + g.filter) # filter out any data outside these bounds s1 = subset(s, (mass_g<(m.fit.raw+g.filter) & mass_g>(m.fit.raw-g.filter))) s1$mass_g.clean = s1$mass_g s1$time_s.clean = s1$time.dec no.masses.clean[row] = length(s1$mass_g.clean) # save the number of masses retained in the final cleaned dataset # check to see that there is enough data to do the fitting # 'degree' must be less than number of unique points; need degree+1 points for fit that is not just the data # don't fit the data if there are no masses before 0815 (495) if ((length(s1$mass_g.clean)<10)|(min(s1$time.dec)>495)) { # Set fit parameters to NA linear.p[row] = NA quadratic.p[row] = NA cubic.p[row] = NA r_sq.fit[row] = NA AIC.fit[row] = NA mass.earliest[row] = NA mass.latest[row] = NA mass.max[row] = NA mass.min[row] = NA mass.dawn[row] = NA mass.7[row] = NA mass.8[row] = NA mass.9[row] = NA mass.10[row] = NA mass.11[row] = NA mass.12[row] = NA mass.13[row] = NA mass.14[row] = NA mass.15[row] = NA mass.16[row] = NA mass.17[row] = NA mass.18[row] = NA mass.19[row] = NA mass.20[row] = NA mass.dusk[row] = NA am.gain[row] = NA time.max.mass[row] = NA time.min.mass[row] = NA time.earliest.mass[row] = NA time.latest.mass[row] = NA row = row + 1 } if ((length(s1$mass_g.clean) >=10)&(min(s1$time.dec)<=495)) { # refit the cleaned dataset m1 = lm(mass_g.clean ~ poly(time_s.clean,degree=3), data = s1, na.action = na.exclude) s1$m.fit = predict(m1) # save the fitted values from the curve fitting # extract p-values from the fit linear.p[row] = summary(m1)$coefficients[2,4] # p-value for the linear term in the orthogonal fit quadratic.p[row] = summary(m1)$coefficients[3,4] # p-value for the quadratic term in the orthogonal fit cubic.p[row] = summary(m1)$coefficients[4,4] # p-value for the cubic term in the orthogonal fit # NB need to fit non-orthogonal polynomial to get the correct coefficients for the fitted line m2 = lm(mass_g.clean ~ time_s.clean + I(time_s.clean^2) + I(time_s.clean^3), data = s1, na.action = na.exclude) # eqn for cubic: y = a + bx + cx^2 + dx^3 where x is the minutes since midnight # extract parameter estimates b, c and d par.a = as.numeric(m2$coefficients[1]) par.b = as.numeric(m2$coefficients[2]) par.c = as.numeric(m2$coefficients[3]) par.d = as.numeric(m2$coefficients[4]) # calculate the fitted mass at each hour by plugging the relevant times into the fitted funtion mass.dawn[row] = par.a + (par.b*375) + (par.c*(375^2)) + (par.d*(375^3)) mass.7[row] = par.a + (par.b*420) + (par.c*(420^2)) + (par.d*(420^3)) mass.8[row] = par.a + (par.b*480) + (par.c*(480^2)) + (par.d*(480^3)) mass.9[row] = par.a + (par.b*540) + (par.c*(540^2)) + (par.d*(540^3)) mass.10[row] = par.a + (par.b*600) + (par.c*(600^2)) + (par.d*(600^3)) mass.11[row] = par.a + (par.b*660) + (par.c*(660^2)) + (par.d*(660^3)) mass.12[row] = par.a + (par.b*720) + (par.c*(720^2)) + (par.d*(720^3)) mass.13[row] = par.a + (par.b*780) + (par.c*(780^2)) + (par.d*(780^3)) mass.14[row] = par.a + (par.b*840) + (par.c*(840^2)) + (par.d*(840^3)) mass.15[row] = par.a + (par.b*900) + (par.c*(900^2)) + (par.d*(900^3)) mass.16[row] = par.a + (par.b*960) + (par.c*(960^2)) + (par.d*(960^3)) mass.17[row] = par.a + (par.b*1020) + (par.c*(1020^2)) + (par.d*(1020^3)) mass.18[row] = par.a + (par.b*1080) + (par.c*(1080^2)) + (par.d*(1080^3)) mass.19[row] = par.a + (par.b*1140) + (par.c*(1140^2)) + (par.d*(1140^3)) mass.20[row] = par.a + (par.b*1200) + (par.c*(1200^2)) + (par.d*(1200^3)) mass.dusk[row] = par.a + (par.b*1095) + (par.c*(1095^2)) + (par.d*(1095^3)) # 'dusk' is 12 hours after dawn mass.dawn.plus2 = par.a + (par.b*495) + (par.c*(495^2)) + (par.d*(495^3)) am.gain[row] = mass.dawn.plus2-mass.dawn[row] # Extract parameters from the fit r_sq.fit[row] = summary(m1)$r.squared # r-squared for the fit AIC.fit[row] = AIC(m1) # AIC for the fit mass.earliest[row] = s1$m.fit[which(s1$time_s.clean==min(s1$time_s.clean))] # the earliest fitted mass mass.latest[row] = s1$m.fit[which(s1$time_s.clean==max(s1$time_s.clean))] # the latest fitted mass mass.max[row] = max(s1$m.fit) # the maximumn fitted mass mass.min[row] = min(s1$m.fit) # the minimum fitted mass time.max.mass[row] = median(s1$time_s.clean[which(s1$m.fit==max(s1$m.fit))]) # finds the time with the max value (median if more than one) time.min.mass[row] = median(s1$time_s.clean[which(s1$m.fit==min(s1$m.fit))]) # finds the time with the min value (median if more than one) time.earliest.mass[row] = min(s1$time_s.clean) # the time at which the first fitted mass occurred time.latest.mass[row] = max(s1$time_s.clean) # the time at which the last fitted mass occurred # we need a rule to exclude any fits when the latest mass that day is substantially before the time of the fit. This is important because # birds sometimes seem to stop foraging part way through the day, especially under FI. # If we implement a rule whereby we exlude any fits if the latest mass of the day was more than 1 hour previously then this appears to elimiate unrealistic fits. # However, we do lose dusk masses for some days. if (time.latest.mass[row] < 480) {mass.9[row]=NA} if (time.latest.mass[row] < 540) {mass.10[row]=NA} if (time.latest.mass[row] < 600) {mass.11[row]=NA} if (time.latest.mass[row] < 660) {mass.12[row]=NA} if (time.latest.mass[row] < 720) {mass.13[row]=NA} if (time.latest.mass[row] < 780) {mass.14[row]=NA} if (time.latest.mass[row] < 840) {mass.15[row]=NA} if (time.latest.mass[row] < 900) {mass.16[row]=NA} if (time.latest.mass[row] < 960) {mass.17[row]=NA} if (time.latest.mass[row] < 1020) {mass.18[row]=NA} if (time.latest.mass[row] < 1080) {mass.19[row]=NA} if (time.latest.mass[row] < 1140) {mass.20[row]=NA} if (time.latest.mass[row] < 1035) {mass.dusk[row]=NA} row = row + 1 s1 = subset(s1, select = c("row.number","time_s.clean","mass_g.clean","m.fit")) d.clean = rbind.data.frame(d.clean, s1) } } } } } d.new = merge(d,d.clean, by = "row.number", all.x = FALSE, all.y = TRUE) fit.results = data.frame(room, bird, treatment, day.exp, day.data, day.treat, no.masses,no.masses.clean, linear.p, quadratic.p, cubic.p, r_sq.fit, AIC.fit, mass.earliest, mass.latest, mass.max, mass.min, mass.dawn, mass.7, mass.8, mass.9, mass.10, mass.11, mass.12, mass.13, mass.14, mass.15, mass.16, mass.17, mass.18, mass.19, mass.20, mass.dusk, am.gain, time.max.mass, time.min.mass, time.earliest.mass, time.latest.mass) write.csv(fit.results, "Exp4_fit_results.csv") describe(fit.results) # clean up the dataframe d.new d.new$mass_g.fit = d.new$m.fit d.new$m.fit = NULL write.csv(d.new, "Exp4_mass_data_fitted_and_cleaned.csv") #### Fit some models to the fitted mass data #### fit.results$treatment = factor(fit.results$treatment, levels = c("FS","FIlow","FIhigh")) # calculate mean baseline dawn mass for each aviary fit.results$aviary = fit.results$room baseline = fit.results[fit.results$treatment=="FS",] baseline.mass = summarise(group_by(baseline, aviary), baseline.mass = mean(mass.dawn, na.rm = TRUE)) baseline.mass$aviary = paste("aviary", baseline.mass$aviary) # dawn mass (mean of fitted mass at 0615) # first create a variable for treatment yesterday fit.results$treatment.yesterday[fit.results$day.treat>1] = as.character(fit.results$treatment[fit.results$day.treat>1]) fit.results$treatment.yesterday = factor(fit.results$treatment.yesterday, levels = c("FS", "FIlow", "FIhigh")) fit.results$treatment.yesterday[fit.results$day.treat==1&fit.results$treatment=="FS"] = "FS" # birds were on FS prior to FS. fit.results$treatment.yesterday[fit.results$day.treat==1&fit.results$treatment=="FIlow"] = "FS" fit.results$treatment.yesterday[fit.results$day.treat==1&fit.results$treatment=="FIhigh"] = "FIlow" # now plot the data fit.results$treatment.yesterday = relevel(fit.results$treatment.yesterday, ref = "FS") mean.fit.mass = summarise(group_by(fit.results,bird,treatment.yesterday), Mean_mass = mean(mass.dawn, na.rm = TRUE)) dawn.mass.pair = summarise(group_by(fit.results,treatment.yesterday,room), mean.mass = mean(mass.dawn, na.rm = TRUE)) #calculate mean mass for each treatment mass.treat = summarise(group_by(mean.fit.mass, treatment.yesterday), mean.mass = mean(Mean_mass), sd.mass = sd(Mean_mass), n.mass = length2(Mean_mass)) # now calculate overall mean for each bird and merge this back into the dataframe mean.mean.mass = summarise(group_by(mean.fit.mass, bird), Mean_mass_bird = mean(Mean_mass)) mean.fit.mass = merge(mean.fit.mass,mean.mean.mass, by="bird") # within-subject centering followed by adding back the grant mean so that phase effects are relocated relative to the grand mean mean.fit.mass$norm.mass = (mean.fit.mass$Mean_mass - mean.fit.mass$Mean_mass_bird)+mean(mean.fit.mass$Mean_mass_bird) fig = ggplot(data=mean.fit.mass, aes(x=treatment.yesterday, y=norm.mass)) + geom_boxplot(aes(fill=treatment.yesterday)) + scale_fill_manual(values=c("cyan", "tomato", "tomato")) + guides(fill=FALSE) + ylim(72.5,82.5) + #ylab("Dawn mass (g)") + ylab("") + xlab("") + scale_x_discrete(labels = c("FS" = "FS", "FIlow"=expression(paste(FI[low])), "FIhigh"=expression(paste(FI[high])))) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + annotate("segment", x=1, xend=2, y=81, yend=81, size=.5) + annotate("text", x=1.5, y=81.25, label = "NS", size=4) + annotate("segment", x=2, xend=3, y=81.5, yend=81.5, size=.5) + annotate("text", x=2.5, y=81.75, label = "NS", size=4) + annotate("segment", x=1, xend=3, y=82, yend=82, size=.5) + annotate("text", x=2, y=82.25, label = "0.1 > p > 0.05", size=4) fig win.metafile("Exp 4 dawn mass.wmf") fig dev.off() m = lmer(mass.dawn ~ treatment.yesterday + (1|room/bird), data=fit.results) summary(m) confint(m) anova(m) # Marginally significant effect of treatment on dawn mass m2 = lmer(scale(mass.dawn) ~ treatment.yesterday + (1|room/bird), data=fit.results) summary(m2) fit.results$treatment.yesterday = relevel(fit.results$treatment.yesterday, ref = "FIlow") m1 = lmer(mass.dawn ~ treatment.yesterday + (1|room/bird), data=fit.results) summary(m1) confint(m1) # output the meta-data output.dawn = data.frame("exp" = "Exp 4","aviary" = "all", "variable" = "mass.dawn", "mean.FS" = mass.treat$mean.mass[1], "mean.FI" = mass.treat$mean.mass[3], "beta" = coef(summary(m))["treatment.yesterdayFIhigh", "Estimate"], "se.beta" = coef(summary(m))["treatment.yesterdayFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2))["treatment.yesterdayFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2))["treatment.yesterdayFIhigh", "Std. Error"]) # Now calculate the effects of FI on dawn mass by aviary fit.results$treatment.yesterday = relevel(fit.results$treatment.yesterday, ref = "FS") # pair 1 fit.results.1 = fit.results[fit.results$room=="1", ] m.1 = lmer(mass.dawn ~ treatment.yesterday + (1|bird), data=fit.results.1) summary(m.1) m2.1 = lmer(scale(mass.dawn) ~ treatment.yesterday + (1|bird), data=fit.results.1) summary(m2.1) output.dawn.1 = data.frame("exp" = "Exp 4","aviary" = "aviary 1", "variable" = "mass.dawn", "mean.FS" = dawn.mass.pair$mean.mass[1], "mean.FI" = dawn.mass.pair$mean.mass[7], "beta" = coef(summary(m.1))["treatment.yesterdayFIhigh", "Estimate"], "se.beta" = coef(summary(m.1))["treatment.yesterdayFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.1))["treatment.yesterdayFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.1))["treatment.yesterdayFIhigh", "Std. Error"]) # pair 2 fit.results.2 = fit.results[fit.results$room=="2", ] m.2 = lmer(mass.dawn ~ treatment.yesterday + (1|bird), data=fit.results.2) summary(m.2) m2.2 = lmer(scale(mass.dawn) ~ treatment.yesterday + (1|bird), data=fit.results.2) summary(m2.2) output.dawn.2 = data.frame("exp" = "Exp 4","aviary" = "aviary 2", "variable" = "mass.dawn", "mean.FS" = dawn.mass.pair$mean.mass[2], "mean.FI" = dawn.mass.pair$mean.mass[8], "beta" = coef(summary(m.2))["treatment.yesterdayFIhigh", "Estimate"], "se.beta" = coef(summary(m.2))["treatment.yesterdayFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.2))["treatment.yesterdayFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.2))["treatment.yesterdayFIhigh", "Std. Error"]) # pair 3 fit.results.3 = fit.results[fit.results$room=="3", ] m.3 = lmer(mass.dawn ~ treatment.yesterday + (1|bird), data=fit.results.3) summary(m.3) m2.3 = lmer(scale(mass.dawn) ~ treatment.yesterday + (1|bird), data=fit.results.3) summary(m2.3) output.dawn.3 = data.frame("exp" = "Exp 4","aviary" = "aviary 3", "variable" = "mass.dawn", "mean.FS" = dawn.mass.pair$mean.mass[3], "mean.FI" = dawn.mass.pair$mean.mass[9], "beta" = coef(summary(m.3))["treatment.yesterdayFIhigh", "Estimate"], "se.beta" = coef(summary(m.3))["treatment.yesterdayFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.3))["treatment.yesterdayFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.3))["treatment.yesterdayFIhigh", "Std. Error"]) # Dusk mass fit.results$treatment = relevel(fit.results$treatment, ref = "FS") # now plot the data splitting FI into two bars # first get the means for each bird in each treatment mean.fit.mass = summarise(group_by(fit.results,bird,treatment), Mean_mass = mean(mass.dusk, na.rm = TRUE)) dusk.mass.pair = summarise(group_by(fit.results,treatment,room), mean.mass = mean(mass.dusk, na.rm = TRUE)) #calculate mean mass for each treatment mass.treat = summarise(group_by(mean.fit.mass, treatment), mean.mass = mean(Mean_mass, na.rm = TRUE), sd.mass = sd(Mean_mass, na.rm = TRUE), n.mass = length2(Mean_mass, na.rm = TRUE)) # now calculate overall mean for each bird and merge this back into the dataframe mean.mean.mass = summarise(group_by(mean.fit.mass, bird), Mean_mass_bird = mean(Mean_mass, na.rm = TRUE)) mean.fit.mass = merge(mean.fit.mass,mean.mean.mass, by="bird") # within-subject centering followed by adding back the grant mean so that phase effects are relocated relative to the grand mean mean.fit.mass$norm.mass = (mean.fit.mass$Mean_mass - mean.fit.mass$Mean_mass_bird)+mean(mean.fit.mass$Mean_mass_bird) fig = ggplot(data=mean.fit.mass, aes(x=treatment, y=norm.mass)) + geom_boxplot(aes(fill=treatment)) + scale_fill_manual(values=c("cyan", "tomato", "tomato")) + guides(fill=FALSE) + #ylab("Dusk mass (g)") + ylab("") + xlab("") + scale_x_discrete(labels = c("FS" = "FS", "FIlow"=expression(paste(FI[low])), "FIhigh"=expression(paste(FI[high])))) + ylim(76,90) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + annotate("segment", x=1, xend=2, y=88, yend=88, size=.5) + annotate("text", x=1.5, y=88.25, label = "0.01 > p > 0.05", size=4) + annotate("segment", x=2, xend=3, y=88.5, yend=88.5, size=.5) + annotate("text", x=2.5, y=88.75, label = "NS", size=4) + annotate("segment", x=1, xend=3, y=89, yend=89, size=.5) + annotate("text", x=2, y=89.25, label = "**", size=4) fig win.metafile("Exp 4 dusk mass.wmf") fig dev.off() # models of dusk mass m = lmer(mass.dusk ~ treatment + (1|room/bird), data=fit.results) summary(m) confint(m) anova(m) # m2 = lmer(scale(mass.dusk) ~ treatment + (1|room/bird), data=fit.results) summary(m2) fit.results$treatment = relevel(fit.results$treatment, ref = "FIlow") m1 = lmer(mass.dusk ~ treatment + (1|room/bird), data=fit.results) summary(m1) confint(m1) # output the meta-data output.dusk = data.frame("exp" = "Exp 4","aviary" = "all", "variable" = "mass.dusk", "mean.FS" = mass.treat$mean.mass[1], "mean.FI" = mass.treat$mean.mass[3], "beta" = coef(summary(m))["treatmentFIhigh", "Estimate"], "se.beta" = coef(summary(m))["treatmentFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2))["treatmentFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2))["treatmentFIhigh", "Std. Error"]) # Now calculate the effects of FI on dusk mass by aviary fit.results$treatment = relevel(fit.results$treatment, ref = "FS") # pair 1 fit.results.1 = fit.results[fit.results$room=="1", ] m.1 = lmer(mass.dusk ~ treatment + (1|bird), data=fit.results.1) summary(m.1) m2.1 = lmer(scale(mass.dusk) ~ treatment + (1|bird), data=fit.results.1) summary(m2.1) output.dusk.1 = data.frame("exp" = "Exp 4","aviary" = "aviary 1", "variable" = "mass.dusk", "mean.FS" = dusk.mass.pair$mean.mass[1], "mean.FI" = dusk.mass.pair$mean.mass[7], "beta" = coef(summary(m.1))["treatmentFIhigh", "Estimate"], "se.beta" = coef(summary(m.1))["treatmentFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.1))["treatmentFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.1))["treatmentFIhigh", "Std. Error"]) # pair 2 fit.results.2 = fit.results[fit.results$room=="2", ] m.2 = lmer(mass.dusk ~ treatment + (1|bird), data=fit.results.2) summary(m.2) m2.2 = lmer(scale(mass.dusk) ~ treatment + (1|bird), data=fit.results.2) summary(m2.2) output.dusk.2 = data.frame("exp" = "Exp 4","aviary" = "aviary 2", "variable" = "mass.dusk", "mean.FS" = dusk.mass.pair$mean.mass[2], "mean.FI" = dusk.mass.pair$mean.mass[8], "beta" = coef(summary(m.2))["treatmentFIhigh", "Estimate"], "se.beta" = coef(summary(m.2))["treatmentFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.2))["treatmentFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.2))["treatmentFIhigh", "Std. Error"]) # pair 3 fit.results.3 = fit.results[fit.results$room=="3", ] m.3 = lmer(mass.dusk ~ treatment + (1|bird), data=fit.results.3) summary(m.3) m2.3 = lmer(scale(mass.dusk) ~ treatment + (1|bird), data=fit.results.3) summary(m2.3) output.dusk.3 = data.frame("exp" = "Exp 4","aviary" = "aviary 3", "variable" = "mass.dusk", "mean.FS" = dusk.mass.pair$mean.mass[3], "mean.FI" = dusk.mass.pair$mean.mass[9], "beta" = coef(summary(m.3))["treatmentFIhigh", "Estimate"], "se.beta" = coef(summary(m.3))["treatmentFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.3))["treatmentFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.3))["treatmentFIhigh", "Std. Error"]) #### Exploring the relationship between body weight and consumption #### # sort f1 by room and day #f1 = f1[order(f1[,1], f1[,4]), ] f1$treatment = relevel(f1$treatment, ref = "FS") f1$room = as.factor(f1$pair) f1 = f1[order(f1$room),] fit.results$day.data.total = rep(seq(1:21),3) mean.bird.mass = summarise(group_by(fit.results, room, treatment.yesterday, day.data.total), bird.dawn.mass = mean(mass.9, na.rm = TRUE), bird.dusk.mass = mean(mass.18, na.rm = TRUE)) mean.bird.mass$row.no = seq(1:length(mean.bird.mass$room)) f1$row.no = seq(1:length(f1$day.exp)) f1 = merge(f1, mean.bird.mass, by = "row.no") f1$eaten.bird = f1$total.eaten/2 # Do treatment and consumption the previous day predict dawn mass? # Note that since food consumption in this experiment is measured at 1100 # Therefore unlike Expts 1-3 we are using the total consumed that day to predict both dawn and dusk mass that day. f1$treatment.yesterday = relevel(f1$treatment.yesterday, ref = "FS") m = lmer(bird.dawn.mass ~ treatment.yesterday + eaten.bird + (1|room.x), data=f1) summary(m) # dawn mass is predicted by consumption. confint(m) anova(m) fig = ggplot(data=f1, aes(x=eaten.bird, y=bird.dawn.mass)) + geom_point(aes(colour=treatment.yesterday)) + geom_smooth(method = "lm", colour = "black") + scale_colour_manual(values=c("cyan", "tomato", "tomato")) + guides(fill=FALSE) + ylab("Dawn mass (g)") + xlab("Food consumed previous day (g)") + #annotate("text", x=2.5, y=85.6, label = "Experiment 3", size=5) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) fig # Do treatment and consumption on the current day predict dusk mass? f1$treatment = relevel(f1$treatment, ref = "FS") m = lmer(bird.dusk.mass ~ treatment + eaten.bird + (1|room.x), data=f1) summary(m) # dusk mass is predicted by treatment that day but not by consumption. confint(m) anova(m) f1$treatment = factor(f1$treatment, levels = c("FS", "FIlow", "FIhigh")) fig = ggplot(data=f1, aes(x=eaten.bird, y=bird.dusk.mass)) + geom_point(aes(colour=treatment)) + geom_smooth(method = "lm", colour = "black") + scale_colour_manual(values=c("cyan", "tomato","tomato")) + guides(fill=FALSE) + ylab("Dusk mass (g)") + xlab("Food consumed that day (g)") + #annotate("text", x=2.5, y=85.6, label = "Experiment 3", size=5) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) fig #### Models of efficiency #### # compute dawn efficiency and ask whether this differs by treatment f1$efficiency.dawn = f1$bird.dawn.mass/f1$eaten.bird mean.treat.room = summarise(group_by(f1, treatment, room.x), mean.efficiency = mean((efficiency.dawn), na.rm=T)) mean.treat = summarise(group_by(mean.treat.room, treatment), mean.efficiency = mean((mean.efficiency), na.rm=T)) # within-based pair centering based on raw values for plotting f1$room = f1$room.x mean.room = summarise(group_by(f1, room), mean.efficiency.room.dawn = mean((efficiency.dawn), na.rm=T)) f1 = merge(f1,mean.room, by="room") # within-pair centering followed by adding back the grand mean so that treatment effects are relocated relative to the grand mean f1$norm.eff = ((f1$efficiency.dawn) - f1$mean.efficiency.room.dawn)+(mean((f1$efficiency.dawn), na.rm=T)) fig = ggplot(data=f1, aes(x=treatment, y=norm.eff)) + geom_boxplot(aes(fill=treatment)) + scale_fill_manual(values=c("cyan", "tomato", "tomato")) + #ylab("Efficiency") + ylab("") + xlab("") + guides (fill = FALSE) + ylim(3,9) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + annotate("segment", x=1, xend=2, y=8, yend=8, size=.5) + annotate("text", x=1.5, y=8.2, label = "NS", size=4) + annotate("segment", x=2, xend=3, y=8.4, yend=8.4, size=.5) + annotate("text", x=2.5, y=8.6, label = "**", size=4) + annotate("segment", x=1, xend=3, y=8.8, yend=8.8, size=.5) + annotate("text", x=2, y=9, label = "***", size=4) fig win.metafile("Exp 4 efficiency.wmf") fig dev.off() # Does treatment predict dawn efficiency? f1$treatment.yesterday = factor(f1$treatment.yesterday, levels = c("FS","FIlow","FIhigh")) f1$treatment.yesterday = relevel(f1$treatment.yesterday, ref = "FS") m = lmer(efficiency.dawn ~ treatment.yesterday + (1|room.x), data=f1, na.action = na.omit) summary(m) confint(m) anova(m) m2 = lmer(scale(efficiency.dawn) ~ treatment.yesterday + (1|room.x), data=f1, na.action = na.omit) summary(m2) f1$treatment.yesterday = relevel(f1$treatment.yesterday, ref = "FIlow") m1 = lmer(efficiency.dawn ~ treatment.yesterday + (1|room.x), data=f1, na.action = na.omit) summary(m1) confint(m1) # output the meta-data output.dawn.eff = data.frame("exp" = "Exp 4","aviary" = "all", "variable" = "efficiency.dawn", "mean.FS" = mean.treat$mean.efficiency[1], "mean.FI" = mean.treat$mean.efficiency[3], "beta" = coef(summary(m))["treatment.yesterdayFIhigh", "Estimate"], "se.beta" = coef(summary(m))["treatment.yesterdayFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2))["treatment.yesterdayFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m))["treatment.yesterdayFIhigh", "Std. Error"]) # Now calculate the effects of FI on dawn efficiency by aviary f1$treatment = relevel(f1$treatment, ref = "FS") # pair 1 f1.1 = f1[f1$room=="1", ] m.1 = lm(efficiency.dawn ~ treatment.yesterday, data=f1.1) summary(m.1) m2.1 = lm(scale(efficiency.dawn) ~ treatment.yesterday, data=f1.1) summary(m2.1) output.dawn.eff.1 = data.frame("exp" = "Exp 4","aviary" = "aviary 1", "variable" = "efficiency.dawn", "mean.FS" = mean.treat.room$mean.efficiency[1], "mean.FI" = mean.treat.room$mean.efficiency[7], "beta" = coef(summary(m.1))["treatment.yesterdayFIhigh", "Estimate"], "se.beta" = coef(summary(m.1))["treatment.yesterdayFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.1))["treatment.yesterdayFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.1))["treatment.yesterdayFIhigh", "Std. Error"]) # pair 2 f1.2 = f1[f1$room=="2", ] m.2 = lm(efficiency.dawn ~ treatment.yesterday, data=f1.2) summary(m.2) m2.2 = lm(scale(efficiency.dawn) ~ treatment.yesterday, data=f1.2) summary(m2.2) output.dawn.eff.2 = data.frame("exp" = "Exp 4","aviary" = "aviary 2", "variable" = "efficiency.dawn", "mean.FS" = mean.treat.room$mean.efficiency[2], "mean.FI" = mean.treat.room$mean.efficiency[8], "beta" = coef(summary(m.2))["treatment.yesterdayFIhigh", "Estimate"], "se.beta" = coef(summary(m.2))["treatment.yesterdayFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.2))["treatment.yesterdayFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.2))["treatment.yesterdayFIhigh", "Std. Error"]) # pair 3 f1.3 = f1[f1$room=="3", ] m.3 = lm(efficiency.dawn ~ treatment.yesterday, data=f1.3) summary(m.3) m2.3 = lm(scale(efficiency.dawn) ~ treatment.yesterday, data=f1.3) summary(m2.3) output.dawn.eff.3 = data.frame("exp" = "Exp 4","aviary" = "aviary 3", "variable" = "efficiency.dawn", "mean.FS" = mean.treat.room$mean.efficiency[3], "mean.FI" = mean.treat.room$mean.efficiency[9], "beta" = coef(summary(m.3))["treatment.yesterdayFIhigh", "Estimate"], "se.beta" = coef(summary(m.3))["treatment.yesterdayFIhigh", "Std. Error"], "standard.beta" = coef(summary(m2.3))["treatment.yesterdayFIhigh", "Estimate"], "standard.se.beta" = coef(summary(m2.3))["treatment.yesterdayFIhigh", "Std. Error"]) #### Final output dataframe for meta-analysis #### output = rbind(output.food, output.dawn, output.dusk, output.dawn.eff) baseline.mean = mean(baseline.mass$baseline.mass) output$baseline.mass = baseline.mean write.csv(output,"Exp 4 data for meta-analysis v3.csv") # output split by aviary output.aviary = rbind(output.eaten.1, output.eaten.2, output.eaten.3, output.dawn.1,output.dawn.2,output.dawn.3, output.dusk.1,output.dusk.2,output.dusk.3, output.dawn.eff.1, output.dawn.eff.2, output.dawn.eff.3) output.aviary = merge(output.aviary,baseline.mass, by = "aviary") write.csv(output.aviary,"Exp 4 data for meta-analysis by aviary v3.csv")