Calculate and plot bootstrapped 95% confidence intervals for model coeficients

# this R code re-produces Figure 1 of … Guerin G.R. & Lowe, A.J. (2013) Leaf morphology shift: new data and analysis support climate link. Biology Letters 9, 20120860. http://rsbl.royalsocietypublishing.org/content/9/1/20120860.full.pdf ... using the data files in the ESM, see http://intl-rsbl.royalsocietypublishing.org/content/9/1/20120860/suppl/DC1

My version of the figure - a) shows coeficients for temporal change in leaf shape (mm/year) using different data splits, and their bootstrapped  95% CIs. b) shows a shifting latitudinal cline in leaf shape between time periods

########begin code
# load data
all <- read.table(file="ESM_angust_all.csv", header=TRUE, sep="\t") #all data
low <- read.table(file="ESM_angust_low_lats.csv", header=TRUE, sep="\t") #data points between -30 and -31 degree lat. post-1950
predvars <- read.table(file="ESM_predvars.csv", header=TRUE, sep="\t") #dummy variables for GLM predictions pre-1950
predvars2 <- read.table(file="ESM_predvars2.csv", header=TRUE, sep="\t") #dummy variables for GLM predictions post-1990

# subset data
mid <- all[9:209,] #1920-1980
post <- all[31:262,] #1950-2011
other <- all[9:262,] #1920-2011
pre <- all[1:30,] #1880-1950
spli <- all[227:262,] #1990-2011
split <- rbind(pre, spli) #combine pre-1950 with post-1990

#set up boots
require(boot) #package must be installed first
bs <- function(formula, data, indices) {
  d <- data[indices,]
  fit <- lm(formula, data=d)
  return(coef(fit)) 
}

bs2 <- function(formula, data, indices) {
  d <- data[indices,]
  fit <- glm(formula, data=d)
  return(coef(fit)) 
}

# run boots and GLM
results <- boot(data=mid, statistic=bs, R=10000, formula=width~Year+latitude)
BOOTY <- boot.ci(results, index=2, type="basic")

results2 <- boot(data=all, statistic=bs, R=10000, formula=width~Year+latitude)
BOOTY2 <- boot.ci(results2, index=2, type="basic")

results3 <- boot(data=post, statistic=bs, R=10000, formula=width~Year+latitude)
BOOTY3 <- boot.ci(results3, index=2, type="basic")

results4 <- boot(data=other, statistic=bs, R=10000, formula=width~Year+latitude)
BOOTY4 <- boot.ci(results4, index=2, type="basic")

results5 <- boot(data=low, statistic=bs, R=10000, formula=width~Year+latitude)
BOOTY5 <- boot.ci(results5, index=2, type="basic")

results6 <- boot(data=pre, statistic=bs, R=10000, formula=width~Year+latitude)
BOOTY6 <- boot.ci(results6, index=2, type="basic")

GLM <- glm(width~latitude + shift, data=split, family=gaussian (link=identity))
pred1 <- predict.glm(GLM, newdata=predvars)
pred2 <- predict.glm(GLM, newdata=predvars2)

results7 <- boot(data=split, statistic=bs2, R=10000, formula=width~shift + latitude)
BOOTY7 <- boot.ci(results7, index=2, type="basic")


# extract results
pre1950 <- c(BOOTY6$t0, BOOTY6$basic[,4:5])

post1950 <- c(BOOTY3$t0, BOOTY3$basic[,4:5])

allx <- c(BOOTY2$t0, BOOTY2$basic[,4:5])

lowx <- c(BOOTY5$t0, BOOTY5$basic[,4:5])

midx <- c(BOOTY$t0, BOOTY$basic[,4:5])

otherx <- c(BOOTY4$t0, BOOTY4$basic[,4:5])

CIs <- as.data.frame(rbind(pre1950, post1950, allx, lowx, midx, otherx))
colnames(CIs) <- c("slope", "lowerCI", "upperCI")
rownames(CIs) <- c("Pre1950", "Post1950", "1880-2011", "North:post1950", "1920-80", "1920-2011")


#plot
nothing <- c(1,2, 3, 4, 5, 6)

par(mar=c(5,4,1,1)+0.1)

layout(matrix(c(1,1,2,2), 2, 2, byrow = TRUE), respect=TRUE, widths=c(2.5,2.5), heights=c(3,3))

plot(CIs$slope~nothing, ylim=c(-0.09, 0.02), cex=1.5, pch=19, ylab="Slope", xlab="", xaxt='n', xlim=c(0.5, 6.5), cex.axis=1, cex.lab=1.5, las=1, bty="l")
arrows(1, CIs[1,1], x1=1, y1=CIs[1, 2], length=0.1, angle=90, lwd=2)
arrows(1, CIs[1,1], x1=1, y1=CIs[1, 3], length=0.1, angle=90, lwd=2)

arrows(2, CIs[2,1], x1=2, y1=CIs[2, 2], length=0.1, angle=90, lwd=2)
arrows(2, CIs[2,1], x1=2, y1=CIs[2, 3], length=0.1, angle=90, lwd=2)

arrows(3, CIs[3,1], x1=3, y1=CIs[3, 2], length=0.1, angle=90, lwd=2)
arrows(3, CIs[3,1], x1=3, y1=CIs[3, 3], length=0.1, angle=90, lwd=2)

arrows(4, CIs[4,1], x1=4, y1=CIs[4, 2], length=0.1, angle=90, lwd=2)
arrows(4, CIs[4,1], x1=4, y1=CIs[4, 3], length=0.1, angle=90, lwd=2)

arrows(5, CIs[5,1], x1=5, y1=CIs[5, 2], length=0.1, angle=90, lwd=2)
arrows(5, CIs[5,1], x1=5, y1=CIs[5, 3], length=0.1, angle=90, lwd=2)

arrows(6, CIs[6,1], x1=6, y1=CIs[6, 2], length=0.1, angle=90, lwd=2)
arrows(6, CIs[6,1], x1=6, y1=CIs[6, 3], length=0.1, angle=90, lwd=2)

abline(h=0, lty=2)

axis(1, at=c(1,2, 3, 4, 5, 6), labels=c("pre","post", "all", "north", "mid", "other"), cex.axis=1.4, las=2)

text(0.55, 0.02, "(a)", cex=2)

plot(width~latitude, data=pre, pch=16, xlab="Latitude", ylab="Leaf width (mm)", cex=1, cex.lab=1.5, cex.axis=1.5, xlim=c(-35, -30), ylim=c(0, 10), bty="l")
new <- all[252:262,]
points(new$latitude, new$width, pch=3, cex=3)
points(spli$latitude, spli$width, pch=1, cex=0.8)

points(predvars$latitude, pred1, type="l", lwd=2, lty=2)
points(predvars2$latitude, pred2, type="l", lwd=2, lty=1)

text(-34.95, 9.8, "(b)", cex=2)

#end