## R script used to analyse the data contained in the datafile lizardFID.csv ## This includes both statistical analyses as well as code used to create the ## graph in the published paper. ## set working directory first dat <- read.csv("lizardFID.csv") summary(dat) ## sample sizes of the different treatments with(dat, table(approach, gaze)) hist(dat$flight_dist_m) ## flight distance rather skewed ## relevel to make approach=tangential the base level dat$approach <- relevel(dat$approach, "tangential") ## relevel to make gaze=averted is already the base level # basic model mod <- lm(flight_dist_m ~ approach*gaze, data=dat) hist(resid(mod)) ## residuals skewed qqnorm(resid(mod)) qqline(resid(mod)) ## residuals highly non-normal ## re-fit model with response variable ln-transformed mod.ln <- lm(log(flight_dist_m) ~ approach*gaze, data=dat) hist(resid(mod.ln)) ## skew of residuals substantialy lessened qqnorm(resid(mod.ln)) qqline(resid(mod.ln)) ## residuals much closer to normal ## first test interaction term in full model anova(mod.ln) ## interaction not significant, main effects appear significant ## estimate main effects and model coeeficients after dropping interaction term mod.ln.add <- update(mod.ln, . ~ . -approach:gaze) anova(mod.ln.add) summary(mod.ln.add) ## coefficients on normal scale: exp(coef(mod.ln.add)) ## draw barplot of means and bootstrapped CIs ## bootstrap for CIs boot.fun <- function(x, N=10000){ out.vec <- numeric(N) for(i in 1:N){ out.vec[i] <- mean(sample(x, replace=TRUE)) } quantile(out.vec, probs=c(0.025, 0.975)) } CIs <- with(dat, tapply(flight_dist_m, list(approach, gaze), boot.fun)) ## means means <- with(dat, tapply(flight_dist_m, list(approach, gaze), mean)) bp <- barplot(means, beside=TRUE, space=c(0.2, 0.8), xlim=c(0.5, 5), ylim=c(0, 5.5), width=0.8, col=c("lightgrey", "white"), ylab="Flight initiation distance (m)", xlab='Gaze', cex.lab=1.7, cex.axis=1.5, cex.names=1.5) abline(h=0,lwd=3) lines(c(bp[1,1], bp[1,1]), CIs[1,1][[1]], lwd=1.75) lines(c(bp[1,2], bp[1,2]), CIs[1,2][[1]], lwd=1.75) lines(c(bp[2,1], bp[2,1]), CIs[2,1][[1]],lwd=1.75) lines(c(bp[2,2], bp[2,2]), CIs[2,2][[1]],lwd=1.75) ## now explore the effect of size-sex category dat$category <- relevel(dat$category, "FYM") mod.cat <- lm(log(flight_dist_m) ~ approach*gaze*category, data=dat) anova(mod.cat) ## drop three-way interaction mod.cat1 <- update(mod.cat, .~. -approach:gaze:category) anova(mod.cat1) ## drop all two-way interactions mod.cat.add <- lm(log(flight_dist_m) ~ approach+gaze+category, data=dat) anova(mod.cat.add) summary(mod.cat.add) ## large lizards have smaller flight distances ## coefficients on normal scale: exp(coef(mod.cat.add)) ## who does the tail raise display? with(dat, table(category, tail_raise)) ## is tail raise by FYM lizards related to approach? with(subset(dat, category=="FYM"), table(approach, tail_raise)) fisher.test(with(subset(dat, category=="FYM"), table(approach, tail_raise))) ## Ans: not strongly ## is tail raise by FYM lizards related to gaze? with(subset(dat, category=="FYM"), table(gaze, tail_raise)) fisher.test(with(subset(dat, category=="FYM"), table(gaze, tail_raise))) ## Ans: somewhat ## does tail raise add additional info on flight distance (for FYM lizards)? mod.tr <- lm(log(flight_dist_m) ~ approach+gaze+tail_raise, data=subset(dat, category=="FYM")) anova(mod.tr) summary(mod.tr) ## exp(coef(mod.tr))