# (d) Competition experiment analysis (including fur removal analyses for sayi and orbicollis when alone) # Does relative distance from carcass vary with treatment, temperature, and relative size of the beetles? #load dataset (d_competition_exp_data.csv) comp.dat <- read.csv(file.choose()) # use dialog box to import data file names(comp.dat) # lists all variable names comp.dat$treatment<-as.factor(comp.dat$treatment) is.factor(comp.dat$treatment) comp.dat$temp<-as.factor(comp.dat$temp) is.factor(comp.dat$temp) hist(comp.dat$ln.dist) hist(comp.dat$ln.mass) comp.dat$block<-as.factor(comp.dat$block) is.factor(comp.dat$block) library(nlme) model.paired<-lme(ln.dist~treatment*temp*ln.mass, random = ~ 1|block, data = comp.dat, method = "REML") # check model plot(model.paired) qqnorm(model.paired) boxplot(resid(model.paired, type="normalized") ~ comp.dat$temp) bartlett.test(resid(model.paired, type="normalized"), comp.dat$temp) boxplot(resid(model.paired, type="normalized") ~ comp.dat$treatment) bartlett.test(resid(model.paired, type="normalized"), comp.dat$treatment) plot(resid(model.paired) ~ comp.dat$ln.mass) hist(resid(model.paired, type="normalized")) shapiro.test(resid(model.paired)) # ok # compare models with different predictors using information-theoretic approach with MuMIn package; rank models by delta AICc library(MuMIn) options(na.action = "na.fail") # prevent fitting models to different comp.datasets model.paired.ML<-lme(ln.dist~treatment*temp*ln.mass, random = ~ 1|block, data = comp.dat, method = "ML") m2<-dredge(model.paired.ML, rank="AICc") subset(m2,delta<2) # best model is saturated anova(model.paired.ML) # 3-way interaction term significant (F = 5.25, p = 0.03) # alone, cool summary(model.paired.ML) # ln.mass slope = -0.186 (SE=1.071) t=-0.17, df=20, P=0.86 # intercept = 0.347 (SE=0.413) t=0.84, df=26, P=0.41 # together, cool comp.dat$treatment2<-NULL comp.dat$treatment2[comp.dat$treatment=="alone"]<-"zalone" comp.dat$treatment2[comp.dat$treatment=="together"]<-"together" model.paired.b<-lme(ln.dist~treatment2*temp*ln.mass, random = ~ 1|block, data = comp.dat, method = "ML") summary(model.paired.b) # ln.mass slope = -0.357 (SE=0.659) t=-0.54, df=20, P=0.59 # intercept = 0.966 (SE=0.387) t=2.49, df=26, P=0.019 (sayi wins more often, regardless of mass) # alone, warm comp.dat$temp2<-NULL comp.dat$temp2[comp.dat$temp=="cool"]<-"zcool" comp.dat$temp2[comp.dat$temp=="warm"]<-"warm" model.paired.c<-lme(ln.dist~treatment*temp2*ln.mass, random = ~ 1|block, data = comp.dat, method = "ML") summary(model.paired.c) # ln.mass slope = -0.154 (SE=0.740) t=-0.21, df=20, P=0.84 # intercept = 0.111 (SE=0.389) t=0.28, df=26, P=0.78 # together, warm model.paired.d<-lme(ln.dist~treatment2*temp2*ln.mass, random = ~ 1|block, data = comp.dat, method = "ML") summary(model.paired.d) # ln.mass slope = -4.343 (SE=0.971) t=-4.471, df=20, P=0.0002 (larger beetle wins, regardless of species) # intercept = 0.240 (SE=0.373) t=0.64, df=26, P=0.53 ############################### # Plot model results library(ggplot2) theme_bw.3<-theme_bw() + theme(legend.position="none") togetherwarmcomp.data <- expand.grid(treatment = "together", temp = "warm", ln.mass = seq(min(comp.dat$ln.mass), max(comp.dat$ln.mass), length = 100), predicted=NA, lwr=NA, upr=NA) mpred.1 <- predict(model.paired.ML, togetherwarmcomp.data, level = 0, se.fit = TRUE) togetherwarmcomp.data$fit <- mpred.1$fit togetherwarmcomp.data$se.fit <- mpred.1$se.fit togetherwarmcomp.data$lwr <- togetherwarmcomp.data$fit - 1.983972*togetherwarmcomp.data$se.fit togetherwarmcomp.data$upr <- togetherwarmcomp.data$fit + 1.983972*togetherwarmcomp.data$se.fit comp.dat.tog.warm<-subset(comp.dat, comp.dat$temp =="warm" & comp.dat$treatment == "together") together.warm <- ggplot(togetherwarmcomp.data, aes(x=ln.mass, y=fit)) + ggtitle("Together, warm") + xlab("Relative size of beetle") + ylab("Relative distance from carcass") + geom_ribbon(aes(ymin=lwr,ymax=upr),fill="gray80",colour=NA) + geom_line(aes(y=fit),size=2,color="red") + geom_abline(intercept = 0, slope = 0, linetype = "dotted") + geom_vline(xintercept = 0, linetype = "dotted") + geom_point(data=comp.dat.tog.warm, shape=1, size=5, aes(x=ln.mass, y=ln.dist, color=treatment)) + scale_color_manual(values="red") + xlim(-1.3,1.3) + ylim(-8.5,8.5) together.warm + theme_bw.3 # figure 2d togethercoolcomp.data <- expand.grid(treatment = "together", temp = "cool", ln.mass = seq(min(comp.dat$ln.mass), max(comp.dat$ln.mass), length = 100), predicted=NA, lwr=NA, upr=NA) mpred.2 <- predict(model.paired.ML, togethercoolcomp.data, level = 0, se.fit = TRUE) togethercoolcomp.data$fit <- mpred.2$fit togethercoolcomp.data$se.fit <- mpred.2$se.fit togethercoolcomp.data$lwr <- togethercoolcomp.data$fit - 1.983972*togethercoolcomp.data$se.fit togethercoolcomp.data$upr <- togethercoolcomp.data$fit + 1.983972*togethercoolcomp.data$se.fit comp.dat.tog.cool<-subset(comp.dat, comp.dat$temp =="cool" & comp.dat$treatment == "together") together.cool <- ggplot(togethercoolcomp.data, aes(x=ln.mass, y=fit)) + ggtitle("Together, cool") + xlab("Relative size of beetle") + ylab("Relative distance from carcass") + geom_ribbon(aes(ymin=lwr,ymax=upr),fill="gray80",colour=NA) + geom_line(aes(y=fit),size=2,color="#00008B") + geom_abline(intercept = 0, slope = 0, linetype = "dotted") + geom_vline(xintercept = 0, linetype = "dotted") + geom_point(data=comp.dat.tog.cool, shape=1, size=5, aes(x=ln.mass, y=ln.dist, color=treatment)) + scale_color_manual(values="#00008B") + xlim(-1.3,1.3) + ylim(-8.5,8.5) together.cool + theme_bw.3 # figure 2c alonewarmcomp.data <- expand.grid(treatment = "alone", temp = "warm", ln.mass = seq(min(comp.dat$ln.mass), max(comp.dat$ln.mass), length = 100), predicted=NA, lwr=NA, upr=NA) mpred.3 <- predict(model.paired.ML, alonewarmcomp.data, level = 0, se.fit = TRUE) alonewarmcomp.data$fit <- mpred.3$fit alonewarmcomp.data$se.fit <- mpred.3$se.fit alonewarmcomp.data$lwr <- alonewarmcomp.data$fit - 1.983972*alonewarmcomp.data$se.fit alonewarmcomp.data$upr <- alonewarmcomp.data$fit + 1.983972*alonewarmcomp.data$se.fit comp.dat.tog.warm<-subset(comp.dat, comp.dat$temp =="warm" & comp.dat$treatment == "alone") alone.warm <- ggplot(alonewarmcomp.data, aes(x=ln.mass, y=fit)) + ggtitle("alone, warm") + xlab("Relative size of beetle") + ylab("Relative distance from carcass") + geom_ribbon(aes(ymin=lwr,ymax=upr),fill="gray80",colour=NA) + geom_line(aes(y=fit),size=2,color="red") + geom_abline(intercept = 0, slope = 0, linetype = "dotted") + geom_vline(xintercept = 0, linetype = "dotted") + geom_point(data=comp.dat.tog.warm, shape=1, size=5, aes(x=ln.mass, y=ln.dist, color=treatment)) + scale_color_manual(values="red") + xlim(-1.3,1.3) + ylim(-8.5,8.5) alone.warm + theme_bw.3 # figure 2b alonecoolcomp.data <- expand.grid(treatment = "alone", temp = "cool", ln.mass = seq(min(comp.dat$ln.mass), max(comp.dat$ln.mass), length = 100), predicted=NA, lwr=NA, upr=NA) mpred.4 <- predict(model.paired.ML, alonecoolcomp.data, level = 0, se.fit = TRUE) alonecoolcomp.data$fit <- mpred.4$fit alonecoolcomp.data$se.fit <- mpred.4$se.fit alonecoolcomp.data$lwr <- alonecoolcomp.data$fit - 1.983972*alonecoolcomp.data$se.fit alonecoolcomp.data$upr <- alonecoolcomp.data$fit + 1.983972*alonecoolcomp.data$se.fit comp.dat.tog.cool<-subset(comp.dat, comp.dat$temp =="cool" & comp.dat$treatment == "alone") alone.cool <- ggplot(alonecoolcomp.data, aes(x=ln.mass, y=fit)) + ggtitle("alone, cool") + xlab("Relative size of beetle") + ylab("Relative distance from carcass") + geom_ribbon(aes(ymin=lwr,ymax=upr),fill="gray80",colour=NA) + geom_line(aes(y=fit),size=2,color="#00008B") + geom_abline(intercept = 0, slope = 0, linetype = "dotted") + geom_vline(xintercept = 0, linetype = "dotted") + geom_point(data=comp.dat.tog.cool, shape=1, size=5, aes(x=ln.mass, y=ln.dist, color=treatment)) + scale_color_manual(values="#00008B") + xlim(-1.3,1.3) + ylim(-8.5,8.5) alone.cool + theme_bw.3 # figure 2a ############################################ # Supplemental analysis with sex combination differences # create a variable with each sex combination for experiments comp.dat$sex.comb <- paste(comp.dat$orb.sex, comp.dat$sayi.sex, sep="_") comp.dat$sex.comb <- as.factor(comp.dat$sex.comb) comp.dat$sex.comb # rerun best-performing model with sex combination model.paired.ML.sex<-lme(ln.dist~treatment*temp*ln.mass + sex.comb*treatment, random = ~ 1|block, data = comp.dat, method = "ML") plot(model.paired.ML.sex) qqnorm(model.paired.ML.sex) boxplot(resid(model.paired.ML.sex, type="normalized") ~ comp.dat$temp) bartlett.test(resid(model.paired.ML.sex, type="normalized"), comp.dat$temp) boxplot(resid(model.paired.ML.sex, type="normalized") ~ comp.dat$treatment) bartlett.test(resid(model.paired.ML.sex, type="normalized"), comp.dat$treatment) boxplot(resid(model.paired.ML.sex, type="normalized") ~ comp.dat$sex.comb) bartlett.test(resid(model.paired.ML.sex, type="normalized"), comp.dat$sex.comb) plot(resid(model.paired.ML.sex) ~ comp.dat$ln.mass) hist(resid(model.paired.ML.sex, type="normalized")) shapiro.test(resid(model.paired.ML.sex)) # ok model.paired.ML<-lme(ln.dist~treatment*temp*ln.mass, random = ~ 1|block, data = comp.dat, method = "ML") AIC(model.paired.ML.sex, model.paired.ML) ############################################ # Fur removed analysis - alone treatment only (sayi and orbicollis) #load dataset (d_fur_removal_data.csv) comp.dat.fur <- read.csv(file.choose()) # use dialog box to import data file names(comp.dat.fur) # lists all variable names comp.dat.fur$species <- as.factor(comp.dat.fur$species) summary(comp.dat.fur$species) comp.dat.fur$sex <- as.factor(comp.dat.fur$sex) summary(comp.dat.fur$sex) summary(comp.dat.fur$temp.C) comp.dat.fur$temp<-NULL comp.dat.fur$temp[comp.dat.fur$temp.C=="5"]<-"cool" comp.dat.fur$temp[comp.dat.fur$temp.C=="15"]<-"warm" comp.dat.fur$temp<-as.factor(comp.dat.fur$temp) is.factor(comp.dat.fur$temp) summary(comp.dat.fur$temp) hist(comp.dat.fur$beetle.mass) hist(comp.dat.fur$mouse.mass) hist(comp.dat.fur$fur.area.removed) comp.dat.fur$ln.fur.removed <- log(comp.dat.fur$fur.area.removed+1) hist(comp.dat.fur$ln.fur.removed) boxplot(comp.dat.fur$ln.fur.removed ~ comp.dat.fur$sex) plot(comp.dat.fur$ln.fur.removed ~ comp.dat.fur$mouse.mass) fur.model <- glm(ln.fur.removed ~ species*temp*beetle.mass, data = comp.dat.fur) plot(fur.model) boxplot(resid(fur.model) ~ comp.dat.fur$temp) bartlett.test(resid(fur.model), comp.dat.fur$temp) boxplot(resid(fur.model) ~ comp.dat.fur$species) bartlett.test(resid(fur.model), comp.dat.fur$species) plot(resid(fur.model) ~ comp.dat.fur$beetle.mass) hist(resid(fur.model)) shapiro.test(resid(fur.model)) # compare models with different predictors using information-theoretic approach with MuMIn package; rank models by delta AICc library(MuMIn) options(na.action = "na.fail") # prevent fitting models to different comp.datasets fur.model <- glm(ln.fur.removed ~ species*temp*beetle.mass, data = comp.dat.fur) mfur<-dredge(fur.model, rank="AICc") subset(mfur,delta<2) # best model has species, temperature, and their interaction best.fur.model <- glm(ln.fur.removed ~ species*temp, data = comp.dat.fur) plot(best.fur.model) boxplot(resid(best.fur.model) ~ comp.dat.fur$temp) bartlett.test(resid(best.fur.model), comp.dat.fur$temp) boxplot(resid(best.fur.model) ~ comp.dat.fur$species) bartlett.test(resid(best.fur.model), comp.dat.fur$species) hist(resid(best.fur.model)) shapiro.test(resid(best.fur.model)) # cool, orbicollis summary(best.fur.model) # at cool temperatures, sayi removed more fur than orbicollis # estimate (sayi vs orbicollis) = 1.098 (SE=0.275), t=3.99, P=0.0002 # orbicollis removed more fur at warm than cool temperatures # estimate (warm vs cool) = 1.817 (SE=0.270), t=6.74, P<0.00001 # cool, sayi comp.dat.fur$species2<-NULL comp.dat.fur$species2[comp.dat.fur$species=="orbicollis"]<-"zorbicollis" comp.dat.fur$species2[comp.dat.fur$species=="sayi"]<-"sayi" best.fur.model.2 <- glm(ln.fur.removed ~ species2*temp, data = comp.dat.fur) summary(best.fur.model.2) # sayi removed more fur at warm than cool temperatures # estimate (warm vs cool) = 1.165 (SE=0.275), t=4.23, P=0.0001 # warm, orbicollis comp.dat.fur$temp2<-NULL comp.dat.fur$temp2[comp.dat.fur$temp=="cool"]<-"zcool" comp.dat.fur$temp2[comp.dat.fur$temp=="warm"]<-"warm" best.fur.model.3 <- glm(ln.fur.removed ~ species*temp2, data = comp.dat.fur) summary(best.fur.model.3) # at warm temperatures, sayi and orbicollis removed similar amounts of fur # estimate (sayi vs orbicollis) = 0.446 (SE=0.270), t=1.66, P=0.10 ################## # Plot raw data library(ggplot2) theme_bw.2<-theme_bw() + theme(legend.direction="horizontal", legend.position=c(0.5,1.047), legend.title=element_blank()) fur.boxplot.raw<-ggplot(comp.dat.fur, aes(y = fur.area.removed, x = species, fill = species)) + facet_grid(~ temp) + geom_boxplot() + scale_fill_manual(values=c("orange1", "mediumpurple4")) fur.boxplot.raw + theme_bw()+theme(legend.position=c(0.20, 0.80)) # figure 3a ################################################################## ##################################################################