--- title: "Skin temperature and reproductive condition in wild female chimpanzees" author: "Guillaume Dezecache, Claudia Wilke, Nathalie Richi, Christof Neumann & Klaus Zuberbühler" date: "19 June 2017" output: pdf_document: toc: true toc_depth: 2 fontsize: 11pt geometry: margin=0.8in --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) options(width=90) ``` ```{r, warning=F, message=FALSE} # load necessary libraries library(lme4); library(car); library(effects); library(MuMIn) sessionInfo() ``` \clearpage # Data set annotation The data for this study are in `chimpskintemperature.csv`. The data set contains the following columns: + **Date** The calendar date the image was recorded. + **fem** Subject ID + **img** Image ID. + **temper** Subject skin temperature ($^\circ$Celsius). + **humid** Humidity (percent). + **dist** Distance between camera and subject (meters). + **ambitemp** Ambient temperature ($^\circ$Celsius). + **bodypart** Bodypart of the subject. + **preg** Was the female pregnant (1) or cycling (0). + **swell** The swelling stage of the female (0 through 4). \clearpage # Read data and raw data figures ```{r, fig.height=7, fig.width=8} xdata <- read.table("chimpskintemperature.csv", header=T, sep=",") mypars <- par(no.readonly = T) # save graphical parameters # settings for graphical adjustments (separation, jitter etc) offs <- 0.15 # separation of groups offs2 <- 0.05 # for point jitter offs3 <- 0.1 # for median width # create layout lmat <- matrix(c(10, 1, 2, 3, 12, 10, 4, 5, 6, 12, 10, 7, 8 ,9, 12, 13, 11, 11, 11, 14), ncol = 5, byrow = T) layout(mat = lmat, widths = c(1, 4, 4, 4, 2), heights = c(4, 4, 4, 1)) # loop through each bodypart for(f in levels(xdata$bodypart)) { # set up empty plot par(mar = c(0.7, 1, 1.7, 0.1)) plot(0, 0, "n", xlim = c(-0.5, 4.5), ylim = c(18, 38), ann = F, axes = F, yaxs = "i", xaxs = "i") # extract relevant data points temp <- xdata[xdata$bodypart == f, ] preg <- temp[temp$preg == 1, ] cycl <- temp[temp$preg == 0, ] # plot points for cycling and swelling x <- preg$swell + offs + runif(n = length(preg$swell), min = offs2*(-1), max = offs2) # add some jitter points(x, preg$temper, pch = 16, col = rgb(0.2, 0.2, 0.2, 0.3), cex = 0.5) x <- cycl$swell - offs + runif(n = length(cycl$swell), min = offs2*(-1), max = offs2) # add some jitter points(x, cycl$temper, pch = 16, col = rgb(1, 0, 0, 0.3), cex = 0.5) # plot medians and quartiles x <- tapply(preg$temper, preg$swell, median) xx <- as.numeric(names(x)) segments(x0 = xx + offs - offs3, y0 = x, x1 = xx + offs + offs3, y1 = x, col = rgb(0, 0, 0), lwd=2) x1 <- tapply(preg$temper, preg$swell, quantile, probs = 0.25) x2 <- tapply(preg$temper, preg$swell, quantile, probs = 0.75) segments(x0 = xx + offs, y0 = x1, x1 = xx + offs, y1 = x2, col = rgb(0, 0, 0), lwd=1) x <- tapply(cycl$temper, cycl$swell, median) xx <- as.numeric(names(x)) segments(x0 = xx - offs - offs3, y0 = x, x1 = xx - offs + offs3, y1 = x, col = rgb(1, 0, 0), lwd = 2) x1 <- tapply(cycl$temper, cycl$swell, quantile, probs = 0.25) x2 <- tapply(cycl$temper, cycl$swell, quantile, probs = 0.75) segments(x0 = xx - offs, y0 = x1, x1 = xx - offs, y1 = x2, col = rgb(1, 0, 0), lwd = 1) # add axes if(f %in% c("Ear", "Foot", "Neck")) { axis(2, tcl = -0.18, xpd = T, las = 1, labels = NA) v <- c(20, 25, 30, 35) text(labels = v, x = -0.68, y = v, xpd = T, cex = 0.8) } if(f %in% c("Neck", "Nose", "Swelling")) { mtext(text = 0:4, side = 1, line = 0.3, xpd = T, at = 0:4, cex = 0.5) } # add body part labels box() rect(-0.5, 38.5, 4.5, 40.5, xpd=T, border = NA, col = "grey") text(2, 39.5, labels = f, xpd = T) } # add overall axes annotations par(mar = c(0, 0, 0, 0)) plot(0, 0, "n", axes = F, ann = F) bq <- bquote(plain("surface temperature (")~degree~plain("Celsius)")) text(-0.1, 0, labels = bq , srt = 90, cex = 1.5) plot(0, 0, "n", axes = F, ann = F) text(0, 0, labels = "swelling stage", cex = 1.5) # legend plot(0,0,"n", axes=F, ann=F) legend("center", col = c("red", "black"), legend = c("cycling", "pregnant"), pch = 16, pt.cex = 1.2, bty = "n") par(mypars) # reset graphical parameters ``` \clearpage # Data preparation ```{r, fig.height=3.2} # check distributions, transform if necessary and standardize tempermin <- min(xdata$temper) # keep minimum temperature (for back transforming in the figures) par(mfcol = c(1, 2)) hist(xdata$temper) hist((xdata$temper - tempermin)^2) xdata$temper <- (xdata$temper - tempermin)^2 hist(xdata$ambitemp) hist((xdata$ambitemp - min(xdata$ambitemp))^2) xdata$ambitemp <- scale(xdata$ambitemp) hist(xdata$dist) hist(sqrt(xdata$dist)) xdata$dist <- scale(sqrt(xdata$dist)) hist(xdata$humid) hist((xdata$humid - min(xdata$humid))^2) xdata$humid <- scale((xdata$humid - min(xdata$humid))^2) # recode swelling as factor xdata$swell <- as.factor(paste0("s", as.character(xdata$swell))) ``` \clearpage # Data analysis ```{r, fig.height=7.5, eval=TRUE} # fit full model res <- lmer(temper ~ swell*preg +ambitemp +dist +humid +(swell +preg|bodypart) +(1|fem/img), data=xdata, REML=F) # diagnostics and assumption checks par(mfrow = c(3, 2)) VarCorr(res) vif(lm(temper ~ swell +preg +ambitemp +dist +humid, data=xdata)) hist(resid(res)) plot(fitted(res), resid(res)) qqnorm(resid(res)); qqline(resid(res)) plot(resid(res) ~ xdata$swell) plot(resid(res) ~ as.factor(xdata$preg)) # fit null model null <- lmer(temper ~ ambitemp +dist +humid +(swell+ preg|bodypart) +(1|fem/img), data=xdata, REML=F) # compare null and full model (likelihood ratio test) anova(null, res, test = "Chisq") # test interaction red <- update(res, .~. -preg:swell) anova(red, res, test = "Chisq") # R2 r.squaredGLMM(res) # show full model results coefficients(summary(res)) ``` \clearpage # Plot model results ```{r} # create plotting data (via effects package) pdatares <- data.frame(effect("swell:preg", res, xlevels = list(preg = c(0, 1)))) # draw empty plot without axes bq <- bquote(plain("surface temperature (")~degree~plain("Celsius)")) plot(0, 0, xlim = c(0.5, 5.5), ylim = c(45, 220), xlab = "swelling stage", ylab = bq, type = "n", axes = F) # make x-coordinates xvals <- c(1:5, 1:5) + rep(c(-0.18, 0.18), each = 5) # add point estimates and CIs points(xvals, pdatares$fit, pch = rep(c(15, 16), each = 5), col = rep(c("red", "black"), each = 5)) arrows(xvals, pdatares$lower, xvals, pdatares$upper, code = 3, angle = 90, length = 0, col = rep(c("red", "black"), each = 5)) box() # draw axes axis(1, at = c(1:5), labels = 0:4) labs <- seq(25, 35, by = 1) ats <- (labs - tempermin)^2 # backtransform to original scale axis(2, las = 1, at = ats, labels = labs) # legend legend("topleft", legend = c("cycling", "pregnant"), col = c( "red", "black"), pch = c(15, 16), horiz = T, bty = "n", pt.cex = 1.3, cex = 0.9) ``` ```{r, fig.height=8, fig.width=8} # create new data frame with predicted values newd <- data.frame(expand.grid(bodypart = levels(xdata$bodypart), swell = levels(xdata$swell), preg = c(0, 1)), humid = 0, dist = 0, ambitemp = 0, fem = NA, img = NA) newd$pred <- predict(res, newd, re.form = ~(swell +preg|bodypart)) # a function to plot per body part plotfunc <- function(bpart, MAIN, YLIM, legendpos = "bottomleft", xdata = newd) { # select bodypart pdata <- xdata[xdata$bodypart == bpart, ] # cycling first, leftish (to the right: pregnant) pdata$xcord <- rep(1:5, 2) + rep(c(-0.15, 0.15), each = 5) # cycling = red, pregnant = black pdata$co <- rep(c("red", "black"), each = 5) # cycling = square, pregnant = circle pdata$ps <- rep(c(15, 16), each = 5) # set up plot and add axes bq <- bquote(plain("surface temperature (")~degree~plain("Celsius)")) plot(0, 0, xlim = c(0.5, 5.5), ylim = YLIM, "n", axes = F, xlab = "swelling stage", ylab = bq, main = MAIN) axis(1, at = 1:5, labels = 0:4,lwd = NA); box() # y-axis: backtransform labs <- seq(25, 35, by = 1) ats <- (labs - tempermin) ^ 2 axis(2, at = ats, labels = labs, las = 1) # separate swelling stages by lines abline(v = seq(0.5, 5.5, by = 1), lty = 3, col = "grey80") # plot point estimates points(pdata$xcord, pdata$pred, pch=pdata$ps, col = pdata$co, cex = 1.5) # add legend legend(legendpos, legend = c("cycling", "pregnant"), pch = c(15, 16), bty = "n", col = c("red", "black"), cex = 0.8) } YLIMS <- c(45, 220) # on the transformed scale par(mfrow = c(3, 3)) plotfunc("Ear", "ear", YLIMS) plotfunc("Eye", "eye", YLIMS) plotfunc("Face", "face", YLIMS) plotfunc("Foot", "foot", YLIMS) plotfunc("Hand", "hand", YLIMS) plotfunc("Mouth", "mouth", YLIMS) plotfunc("Neck", "neck", YLIMS) plotfunc("Nose", "nose", YLIMS) plotfunc("Swelling", "swelling", YLIMS) ```