library(smatr) library(lattice) library(nlme) library(ade4) p.rounder <- function(x, alpha = 0.05, stars = T){ y <- numeric(length(x)) for(i in 1:length(x)) { if(stars != T){ y[i] = "ns" if(x[i] < alpha) {y[i] = sprintf("%.4f", x[i])} if(x[i] < 0.0001){y[i] = "< 0.0001"} } else { y[i] = "ns" if(x[i] < 0.05) {y[i] = "*"} if(x[i] < 0.001) {y[i] = "**"} if(x[i] < 0.0001){y[i] = "***"} } } return(y) } # Read in data ind <- read.csv("data/ind_2009-08-18.csv", sep = ";", skip = 2) spp <- read.csv("data/tax_2009-10-23.csv", sep = ";", skip = 2) sap <- read.csv("data/Physiology/Sapflow ALL COMBINED.csv") def <- read.csv("data/Herbivory/bark_defense 2008 08 27.csv") roth<- read.csv("data/anatomy/table 1-48 of Roth 1981.csv") phy <- read.tree("data/phylo 2009 09 16.txt") stiff<-read.csv("data/stability/Modulus-Bark-Wood 2010 01 29.csv") # clean up data a bit & make a few necessary variables # IND ind$fire <- spp$fire[match(ind$taxon_code, spp$taxon_code)] ind$ba <- pi * (ind$DBH/200)^2 ind <- ind[ind$bark_thick < 40,] # drop out the super-thick values of Manilkara COLS <- c("Family", "Genus", "species", "bar_code", "plot_code", "DBH", "bark_thick", "twig_diam", "twig_bark_thick", "twig_dens", "sapwood_dens", "thickness", "fire") dat <- ind[!is.na(ind$bark_thick), which(names(ind) %in% COLS)] # slim dataset down to relevant columns, and only indivs with bark measurements dat$twig_bark_thick[dat$bar_code == "T10424"] <- NA # Eliminate a bad twig barkmeasurement dat <- dat[dat$Family != "Arecaceae",] # drop out all palms - no bark there dat$trunk_perc <- 100 * dat$bark_thick/dat$DBH/2/10 # calculate bark thickness as a fraction of trunk diam dat$twig_perc <- 100 * dat$twig_bark_thick/dat$twig_diam/2 # calculate bark thickness as a fraction of twig diam dat$Genus <- as.character(dat$Genus) dat$Genus[dat$Genus == "IND"] <- paste(as.character(dat$Family[dat$Genus == "IND"]), "IND", sep = "") dat$Genus <- factor(dat$Genus) dat$name <- paste(dat$Genus, dat$species) # OTHER non_endemic <- dat[dat$fire>=0.14,] # seperate out the species that are only found in low-fire habitats plots <- sort(unique(ind$plot_code)) fams <- table(dat$Family, exclude = "Arecaceae") fams <- -sort(-fams) #SPP spp$name <- paste(spp$Genus, spp$species) # STIFF stiff.trim <- data.frame( Family = stiff$Family, Genus = stiff$Genus, species = stiff$Species, name = stiff$name, ind = stiff$ind, form = stiff$SimpleLifeForm, part = stiff$SimpleSegment, wood.thick = stiff$Diameter.W/2, wood.modulus = stiff$Module.W, bark.modulus = stiff$E.B, bark.thick = (stiff$Diameter.WB - stiff$Diameter.W)/2, perc.stiff = (stiff$EI.WB - stiff$EI.W)/stiff$EI.W ) stiff.trim$perc.thick <- stiff.trim$bark.thick/stiff.trim$wood.thick stiff.trim$perc.modulus <- stiff.trim$bark.modulus/stiff.trim$wood.modulus stiff.trim$COL <- as.numeric(stiff.trim$form) stiff.trim$PCH <- as.numeric(stiff.trim$part) stiff.trim <- stiff.trim[stiff.trim$form == "Tree",] stiff.mean <- aggregate(stiff.trim[,9:14], by = list(part = stiff.trim$part, ind = stiff.trim$ind), mean, na.rm = T) stiff.sd <- aggregate(stiff.trim[,9:14], by = list(part = stiff.trim$part, ind = stiff.trim$ind), sd) stiff.n <- aggregate(stiff.trim[,9:14], by = list(part = stiff.trim$part, ind = stiff.trim$ind), length) stiff.se <- cbind(stiff.mean[,1:2], stiff.sd[,3:8]/sqrt(stiff.n[,3:8])) stiff.se[is.na(stiff.se)] <- 0 stiff.mean$PCH <- as.numeric(stiff.mean$part) # DEF temp <- aggregate(dat$bark_thick, by = list(name = dat$name),mean, na.rm = T) def$bark_thick <- temp$x[match(def$name, temp$name)] def$quantity <- factor(def$quantity, ordered = T) ################ ### TABLE 1 ### ################ # Fnctonal hypotheses - made in Word. ################ ### TABLE 2 ### ################ # Partition variance in bark thickness among taxonomic levels tax.par1 <- lme(bark_thick ~ 1, random = ~1|Family/Genus/species, data = dat, na.action = na.exclude, subset = DBH > 20) summary(tax.par1) vars <- c(fam = 1.387883, gen = 1.148435, spp = 1.713495, ind = 3.287942)^2 round(100* vars/sum(vars), 2) #same for twigs tax.par.twig <- lme(twig_bark_thick ~ 1, random = ~1|Family/Genus/species, data = dat, na.action = na.exclude) summary(tax.par.twig) vars <- c(fam = 0.4745464, gen = 0.2458166, spp = 0.5376021, ind = 0.8304573)^2 round(100* vars/sum(vars), 2) ################ ### FIGURE 1 ### ################ # Panel A # trunk bark v DBH in overall pdf(file = "figs/Figure 1 2010 05 13.pdf", paper = "usr", width = 11, height = 8.5) par(mar = c(4, 4, 1, 0), oma =c(0, 0, 0, 0), las = 1, bty = "n", pty = "s", tcl = 0.5) YLIM = c(0.5, 40) layout(mat = matrix(c(1, 2, 3, 3, 4, 5), nrow = 2), widths = c(1.0, 0.9, 1.0)) Xes <- data.frame(DBH = seq(10, max(dat$DBH, na.rm = T), by = 2)) plot(dat$DBH, jitter(dat$bark_thick, amount = 0.1), xlab = "DBH (cm)", ylab = "Trunk bark thickness (mm)", ylim = YLIM, xlim = c(10, 175), axes = F, pch = 16, col = "gray", cex = 0.6) axis(1, at = c(10, 50, 100, 150)) axis(2) lm1 <- lm(bark_thick ~ DBH, data = dat) dat$resid[!is.na(dat$bark_thick) & !is.na(dat$DBH)] <- resid(lm1) Asym <- 10 K <- 30 temp <-na.omit(data.frame(bark_thick = dat$bark_thick, DBH = dat$DBH)) asy2 <- nls(bark_thick ~ SSmicmen(DBH, Asym, K), data = temp) P.lm <- bquote(paste("Linear: P < 0.0001, ", AIC == .(sprintf("%1.0f", AIC(lm1))))) P.asy <- bquote(paste("Asymptotic: P < 0.0001, ", AIC == .(sprintf("%1.0f", AIC(asy2))))) pred.lm <- predict(lm1, newdata = Xes) pred.asy <- predict(asy2, newdata = Xes) lines(Xes$DBH, pred.lm) lines(Xes$DBH, pred.asy, col = "blue", lty = 2) abline(h = coef(asy2)[1], col = "blue", lty = 3) text(10, 37, P.lm, pos = 4, family = "mono") text(10, 35, P.asy, pos = 4, family = "mono", col = "blue") text(10, 40, "a)", pos = 4, family = "Times", cex = 1.3) mtext("Figure 1", adj = 0.05, family = "Times", cex = 1.3, outer = T, line = -2) # Panel B # twig bark v twig diam overall YLIM = c(0.01, 8) Xes <- data.frame(twig_diam = seq(min(dat$twig_diam, na.rm = T), max(dat$twig_diam, na.rm = T), by = 0.1)) plot(dat$twig_diam, jitter(dat$twig_bark_thick, amount = 0.1), xlab = "Twig diameter (mm)", ylab = "Twig bark thickness (mm)", ylim = YLIM, xlim = c(2, 27), pch = 16, col = "gray", cex = 0.6) lm1 <- lm(twig_bark_thick ~ twig_diam, data = dat) lm2 <- lm(log(twig_bark_thick) ~ twig_diam, data = dat) P.lm <- bquote(paste("Linear: P < 0.0001, ", AIC == .(sprintf("%1.0f", AIC(lm1))))) P.lm2 <- bquote(paste("Logarithmic: P < 0.0001, ", AIC == .(sprintf("%1.0f", AIC(lm2))))) lines(Xes$twig_diam, predict(lm1, newdata = Xes), col = "black") lines(Xes$twig_diam, exp(predict(lm2, newdata = Xes)), col = "red", lty = 2) text(2, 7.5, P.lm, pos = 4, family = "mono") text(2, 7.1, P.lm2, pos = 4, family = "mono", col = "red") text(2, 8, "b)", pos = 4, family = "Times", cex = 1.3) # Panel C #Bark thickness in all families with at least 10 individuals par(mar = c(4, 8, 1, 1), las = 1, bty = "n", pty = "m") temp <- dat[dat$Family %in% names(fams)[fams > 10] & dat$DBH < 20,] temp$Family <- factor(temp$Family) bymedian <- reorder(temp$Family, temp$trunk_perc, median, na.rm = T) plot.new() plot.window(xlim = c(0.1, 7), ylim = c(1, 36), log = "x") boxplot(temp$trunk_perc ~ bymedian, xlab = "Bark thickness (% of trunk radius)", horizontal = T, add = T, col = "gray", boxlty = 0, whisklty = 1, medlwd = 1, axes = F) axis(1, at = c(0.1, 0.2, 1, 5), labels = c(0.1, 0.2, 1, 5)) axis(2, at = 1:36, labels = names(sort(attr(bymedian, "scores"))), lwd = 0) par(xpd = NA) #text(0.05, 1:36, names(sort(attr(bymedian, "scores"))), pos = 2) segments(0.01, seq(1.5, 35.5), 5, seq(1.5, 35.5), col = "gray80") text(0.01, 37, "c)", family = "Times", cex = 1.3) par(xpd = F) # Panel D: Trunk % par(mar = c(4, 4, 1, 0)) w = 1/6 hist.all <- hist(dat$trunk_perc, breaks = seq(0, 7, 2*w), plot = F) hist.end <- hist(non_endemic$trunk_perc, plot = F, breaks = seq(0, 7, 2*w)) # ks.test(x = hist.all$density[1:19], y = hist.end$density[1:19]) plot.new() plot.window(xlim = c(0, 5), ylim = c(0, 0.4)) axis(1); axis(2) title(xlab = "Bark thickness (% of trunk radius)", ylab = "Frequency") rect(hist.all$mids - w, 0, hist.all$mids + w, hist.all$counts/sum(hist.all$counts), col = "gray", border = "transparent") rect(hist.end$mids - w, 0, hist.end$mids + w, hist.end$counts/sum(hist.all$counts), col = "gray20", border = "transparent") legend("topright", legend = c("All species", "Species that range\ninto fire-prone habitats"), fill = c("gray", "gray20"), bty = "n", inset = 0.1, cex = 1.1, border = NA) text(0, 0.4, "d)", family = "Times", pos = 4, cex = 1.3) # Panel E: Twig % w = 1 hist.all <- hist(dat$twig_perc, breaks = seq(0, 35, 2*w), plot = F) hist.end <- hist(non_endemic$twig_perc, plot = F, breaks = seq(0, 35, 2*w)) # ks.test(x = hist.all$density[1:13], y = hist.end$density[1:13]) plot.new() plot.window(xlim = c(0, 30), ylim = c(0, 0.25)) axis(1); axis(2) title(xlab = "Bark thickness (% of twig radius)", ylab = "Frequency") rect(hist.all$mids - w, 0, hist.all$mids + w, hist.all$counts/sum(hist.all$counts), col = "gray", border = "transparent") rect(hist.end$mids - w, 0, hist.end$mids + w, hist.end$counts/sum(hist.all$counts), col = "gray20", border = "transparent") text(0, 0.25, "e)", family = "Times", pos = 4, cex = 1.3) dev.off() # summary stats for big BRIDGE dataset... median(dat$twig_bark_thick, na.rm = T) # 1.917 mean(dat$twig_bark_thick, na.rm = T) # 2.078 max(dat$twig_bark_thick, na.rm = T) # 7.307 median(dat$bark_thick, na.rm = T) # 4 mean(dat$bark_thick, na.rm = T) # 4.53 max(dat$bark_thick, na.rm = T) # 29 mean (dat$bark_thick[dat$DBH >= 20], na.rm = T) # 6.37 range(dat$bark_thick[dat$DBH >= 20], na.rm = T) # 0.5-29 median(dat$twig_perc, na.rm = T) # 8.54 mean(dat$twig_perc, na.rm = T) # 8.97 max(dat$twig_perc, na.rm = T) # 25.46 median(dat$trunk_perc, na.rm = T) # 0.92 mean(dat$trunk_perc, na.rm = T) # 1.01 max(dat$trunk_perc, na.rm = T) # 6.15 pout <- dat[dat$Genus == "Pouteria" & dat$DBH < 20,] hist(pout$bark_thick) range(pout$bark_thick, na.rm = T) inga <- dat[dat$Genus == "Inga" & dat$DBH < 20,] hist(inga$bark_thick) range(inga$bark_thick) aggregate(dat$trunk_perc[dat$DBH < 20], by = list(dat$Family[dat$DBH < 20]), median, na.rm = T) # additional test (reported in text) # does twig bark thickness (as fraction of diameter) change with diameter? lm.twig_diam <- lm(twig_perc ~ twig_diam, data = dat) summary(lm.twig_diam) range(dat$twig_diam) #plot(dat$twig_diam, dat$twig_perc) #abline(lm.twig_diam) ################ ### FIGURE 2 ### ################ # Panel A: Association with Fire-prone habitats pdf(file = "figs/Figure 2 2010 05 13.pdf", paper = "us", width = 8.5, height = 11) par(mfrow = c(1, 2), las = 1, bty = "n", oma = c(1, 3, 1, 1), pty = "s", mar = c(3, 6, 0, 0), tcl = 0.5, pin = c(3, 3), xpd = NA) temp <- aggregate(dat$resid, by = list(name = dat$name), mean, na.rm = T) spp$resid <- temp$x[match(spp$name, temp$name)] #all spp lm.fire2 <- lm(resid ~ log(fire), data = spp, na.action = na.exclude) AIC(lm.fire2) anova(lm.fire2) summary(lm.fire2) YLIM = c(-10, 20) new <- data.frame(fire = seq(0.013, 7.29, 0.1)) pred <- predict(lm.fire2, new, se.fit = T) P <- bquote(paste("All species: ", P== .(sprintf("%1.2f", summary(lm.fire2)$coef[2, 4])))) plot(spp$fire, spp$resid, log = "x", ylim = YLIM, ylab = "Residuals from thickness-DBH relationship (mm)", xlab = "Fire-association index", axes = F, pch = 16, col = "gray", xlim = c(0.01, 7.29)) axis(1, at = c(0.01, 0.1, 1, 5), labels = c(0.01, 0.1, 1, 5)) axis(2) #, labels = NA) lines(new$fire, pred$fit) lines(new$fire, pred$fit + pred$se.fit, lty = 2) lines(new$fire, pred$fit - pred$se.fit, lty = 2) text(0.01, 20, "a)", family = "Times", pos = 4, cex = 1.3) text(0.01, 17, P, pos = 4, cex = 1) # fire-assocaited spp temp <- spp[spp$fire > 0.014,] lm.fire3 <- lm(resid ~ log(fire), data = temp, na.action = na.exclude) summary(lm.fire3) new <- data.frame( fire = seq(min(temp$fire, na.rm = T), max(temp$fire, na.rm = T), length = 50)) pred <- predict(lm.fire3, new, se.fit = T) P <- bquote(paste("Species associated\nwith fire-prone habitats: ", P== .(sprintf("%1.3f", summary(lm.fire3)$coef[2, 4])))) lines(new$fire, pred$fit, col = "red") lines(new$fire, pred$fit + pred$se.fit, lty = 2, col = "red") lines(new$fire, pred$fit - pred$se.fit, lty = 2, col = "red") text(0.01, 12, P, pos = 4, cex = 1, col = "red") mtext("Figure 2", adj = 0, family = "Times", cex = 1.3, line = 2) # HERBIVORE DEFENSE # Panel B: Herbivore defense - Bark thickness versus latex quantity # lok for tradeoffs interspecific between latex quantity, VOC richness/diversity and bark thickness # first, estimate bark thickness for a tree of DBH=10 for all species in the mono dataset #par(mar = c(3, 2, 0, 0)) cor(def[,c(5, 8, 10)], use = "pair") levels(def$quantity) <- c("0", "1", "2", "3", "3") def$resid <- spp$resid[match(def$name, spp$name)] lm.def <- aov(resid ~ quantity * mono_rich, data = def) lm.def1 <- aov(resid ~ quantity + mono_rich, data = def) #AIC(lm.def) #AIC(lm.def1) #summary(lm.def1) #anova(lm.def1) X <- def$quantity Y <- def$resid P <- bquote(paste(P== .(sprintf("%1.4f", summary(lm.def1)[[1]][1,5])))) plot.new() plot.window(xlim = c(0.5, 4.5), ylim = YLIM) title(xlab = "Latex quantity", line = 5) #plot(X, Y, pch = 16, type = "n", col = "gray", axes = F, ylim = YLIM, xlab = "Latex quantity", ylab = "") # plot the points axis(1, at = 1:4, labels = c("None", "Scant", "Moderate", "Abundant"), las = 2) axis(2, labels = NA) new <- data.frame(quantity = factor(seq(0, 3, 1)), mono_rich = mean(def$mono_rich, na.rm = T)) pred <- predict(lm.def1, new, se.fit = T) points(jitter(1 + as.numeric(as.character(def$quantity)), amount = 0.1), Y, pch = 16, col = "gray") points(new$quantity, pred$fit, pch = 16, cex = 1.3) arrows(1:4, pred$fit - pred$se.fit, 1:4, pred$fit + pred$se.fit, lty = 1, angle = 90, length = 0.05, code = 3) text(1, 17, P, pos = 4, cex = 1) text(1, 20, "b)", family = "Times", pos = 4, cex = 1.3) dev.off() ################ ### FIGURE 3 ### ################ # Niklas (1999) says that if bark stiffness is 1/2 of wood stiffness, then bark 32% as thick as the wood radius will contribute equally to overall branch stiffness # mean twig bark thickkness, as a percent of twig radius is 8.97% # According to Niklas's model, twig bark stiffness would need to be >> than that of wood for each to have an equal contribution to twig stiffness # like Niklas's Figure 7a, but with our data, and transposed pred <- function(x, multiple){(1+multiple/x)^0.25-1} pdf(file = "figs/Figure 3 STIFFNESS 2010 04 10.pdf", paper = "us", width = 8.5, height = 11) par(mfrow = c(1, 1), las = 1, bty = "n", pty = "m", pch = 16, pin = c(3, 3), tcl = 0.5) plot.new() plot.window(xlim = c(-0.12, 0.5), ylim = c(0, 0.5)) title(xlab = "Bark stiffness/wood stiffness", ylab = "Bark thickness/wood thickness") Xes <- seq(0, .5, 0.001) lines(Xes, pred(Xes, 1/100), col = "gray") lines(Xes, pred(Xes, 1/10), col = "gray") lines(Xes, pred(Xes, 1/5), col = "gray") lines(Xes, pred(Xes, 1/2), col = "gray") lines(Xes, pred(Xes, 1/1), col = "gray") par(xpd = T) text(x = 0.5, y = pred(0.5, c(1/100, 1/10, 1/5, 1/2, 1)), c("1:100", "1:10", "1:5", "1:2", "1:1"), cex = 0.8, pos = 4) par(xpd = F) with(stiff.mean, points(perc.modulus, perc.thick, col = c("gray", "black")[PCH])) axis(2, at = seq(0, 0.5, 0.1)) axis(1, at = seq(0, 0.5, 0.1)) legend(x = 0.25, y = 0.48, legend = c("Branch", "Trunk"), pch = 16, col = c("gray", "black"), bg = "white") mtext("Figure 3", adj = 0.05, family = "Times", cex = 1.3, line = 2) dev.off() mean(stiff.trim$bark.modulus/stiff.trim$wood.modulus) # 0.1211090 mean(stiff.trim$bark.thick/stiff.trim$wood.thick) # 0.1851890 # Mean contribution of bark to flexural rigidity: stiff.mean$contrib.bark <- stiff.mean$perc.modulus *((stiff.mean$perc.thick + 1)^4-1) mean(stiff.mean$contrib.bark) # 0.1078825 mean(stiff.mean$contrib.bark[stiff.mean$part == "Branch"]) # 0.1130841 mean(stiff.mean$contrib.bark[stiff.mean$part == "Trunk"]) # 0.09979112 ##### Discussion ################ ### FIGURE 4 ### ################ # Analyses of Roth 1981 table # N. periderms predicts depth of rhitdome. and depth of rhitidome predicts total bark thickness lm.bark <- lm(entire.bark.width ~ DBH, dat = roth) summary(lm.bark) roth$resid[!is.na(roth$entire.bark.width) & !is.na(roth$DBH)] <- resid(lm.bark) mean(roth$entire.bark.width[roth$DBH > 20], na.rm = T) sd(roth$entire.bark.width[roth$DBH > 20], na.rm = T)/sqrt(nrow(roth[roth$DBH > 20,])) mean(dat$bark_thick[dat$DBH > 20], na.rm = T) sd(dat$bark_thick[dat$DBH > 20], na.rm = T)/sqrt(nrow(dat[dat$DBH > 20,])) pdf(file = "figs/Figure 4 2010 05 13.pdf", paper = "us", width = 8.5, height = 11) par(mfrow = c(1, 3), las = 1, bty = "n", pty = "s", tcl = 0.5) YLIM = c(-20, 30) XLIM = c(0.8, 25) plot(roth$rhytidome.width+1 , roth$resid, ylab = "Residuals from thickness-DBH relationsip (mm)", xlab = "Rhytidome thickness (mm)", pch = 16, col = "gray", log = "x", ylim = YLIM, xlim = XLIM) lm.peri <- lm(resid ~ log(rhytidome.width+1), data = roth) Xes <- data.frame(rhytidome.width = 1+seq(min(roth$rhytidome.width, na.rm = T), max(roth$rhytidome.width, na.rm = T), 0.1)) P <- bquote(paste("P < 0.0001, ", R^2== .(sprintf("%1.2f", summary(lm.peri)$r.squared)))) pred <- predict(lm.peri, newdata = Xes, se.fit = T) lines(Xes$rhytidome.width, pred$fit) lines(Xes$rhytidome.width, pred$fit + pred$se.fit, lty = 2) lines(Xes$rhytidome.width, pred$fit - pred$se.fit, lty = 2) text(1, 25, P, pos = 4) text(1, 30, "a)", family = "Times", pos = 4, cex = 1.3) mtext("Figure 4", adj = 0, family = "Times", cex = 1.3, line = 2) plot(roth$phelloderm.width+0.01, roth$resid, xlab = "Phelloderm thickness (mm)", ylab = "Residuals from thickness-DBH relationsip (mm)", pch = 16, col = "gray", log = "x", ylim = YLIM, axes = F) axis(1, at = c(0.01, 0.1, 1, 10), labels =c(0.01, 0.1, 1, 10)) axis(2) lm.peri <- lm(roth$resid ~ log(phelloderm.width+0.01), data = roth) Xes <- data.frame(phelloderm.width = 0.01 + seq(min(roth$phelloderm.width, na.rm = T), max(roth$phelloderm.width, na.rm = T), 0.1), DBH = 10) P <- bquote(paste(P== .(sprintf("%1.2f", summary(lm.peri)$coef[2, 4])))) pred <- predict(lm.peri, newdata = Xes, se.fit = T) lines(Xes$phelloderm.width, pred$fit) lines(Xes$phelloderm.width, pred$fit + pred$se.fit, lty = 2) lines(Xes$phelloderm.width, pred$fit - pred$se.fit, lty = 2) text(0.01, 25, P, pos = 4) text(0.01, 30, "b)", family = "Times", pos = 4, cex = 1.3) plot(roth$n.periderm.layers, roth$rhytidome.width+1, ylab = "Rhytidome thickness (mm)", xlab = "Number of phyllogens", pch = 16, col = "gray", log = "xy", xlim = XLIM, ylim = XLIM) lm.peri <- lm(log(roth$rhytidome.width+1) ~ log(n.periderm.layers) + DBH, data = roth) Xes <- data.frame(n.periderm.layers = seq(min(roth$n.periderm.layers, na.rm = T), max(roth$n.periderm.layers, na.rm = T), 0.1), DBH = 10) P <- bquote(paste("P < 0.0001, ", R^2== .(sprintf("%1.2f", summary(lm.peri)$r.squared)))) pred <- predict(lm.peri, newdata = Xes, se.fit = T) lines(Xes$n.periderm.layers, exp(pred$fit)) lines(Xes$n.periderm.layers, exp(pred$fit) + exp(pred$se.fit), lty = 2) lines(Xes$n.periderm.layers, exp(pred$fit) - exp(pred$se.fit), lty = 2) text(1, 20, P, pos = 4) text(1, 25, "c)", family = "Times", pos = 4, cex = 1.3) dev.off() # SUPPLEMENTAL FIGURE 1 # trunk bark v DBH in top 30 spp spp <- table(dat$name) spp <- -sort(-spp) XLIM = c(10, 120) YLIM = range(0, 40) pdf(file = "figs/Figure Supplemental 1 2010 01 31.pdf", paper = "usr", width = 11, height = 8.5) par(mfrow = c(6, 5), mar = c(1, 1, 1, 0), pty = "m", oma =c(3, 3, 0, 0), las = 1, tcl = 0.5) for(i in 1:30){ dat.i <- dat[dat$name == names(spp)[i],] Xes.i <- data.frame(DBH = seq(10, max(dat.i$DBH, na.rm = T), by = 1)) plot(dat.i$DBH, dat.i$bark_thick, ylab = "", xlab = "", xlim = XLIM, ylim = YLIM, axes = F, col = "gray", pch = 16) if(i %% 5 == 1){axis(2)} else{axis(2, labels = NA)} if(i > 25) {axis(1)} else{axis(1, labels = NA)} lm1 <- lm(bark_thick ~ DBH, data = dat.i, na.action = na.exclude) print(summary(lm1)$coef[1,1:2]) stat.lin <- bquote(paste("Lin: ", R^2== .(sprintf("%1.2f", cor(dat.i$bark_thick, fitted(lm1))^2)), " ", AIC == .(sprintf("%1.0f", AIC(lm1))))) lines(Xes.i$DBH, predict(lm1, newdata = Xes.i)) text(10, 39, paste(names(spp)[i], " (", nrow(dat.i), ")", sep = ""), font = 3, pos = 4) text(10, 35, stat.lin, pos = 4, family = "mono") Asym <- 10 K <- 50 temp <-na.omit(data.frame(bark_thick = dat.i$bark_thick, DBH = dat.i$DBH)) asy1 <- NA try(asy1 <- nls(bark_thick ~ SSmicmen(DBH, Asym, K), data = temp, na.action = na.omit), silent = T) if(!is.na(asy1)){ stat.asy <- bquote(paste("Asy.: ", R^2== .(sprintf("%1.2f", cor(dat.i$bark_thick, fitted(asy1))^2)), " ", AIC == .(sprintf("%1.0f", AIC(asy1))))) pred.asy <- predict(asy1, newdata = Xes.i) abline(h = coef(asy1)[1], col = "blue", lty = 2) text(10, 31, stat.asy, pos = 4, family = "mono", col = "blue") lines(Xes.i$DBH, pred.asy, col = "blue") } } title(xlab = "DBH (cm)", ylab = "Trunk bark thickness (mm)", outer = T, line = 1.5, cex.lab = 1.3) mtext("Supplemental Figure 1", adj = 0.05, family = "Times", cex = 1.1, outer = T, line = -1) dev.off() # SUPPLEMENTAL FIGURE 2 # Herbvore defense - Bark thickness versus diversity of monoterpenes pdf(file = "figs/Figure Supplemental 2 2010 05 13.pdf", paper = "us", width = 11, height = 8.5) par(mfrow = c(1, 1), las = 1, bty = "n", oma = c(1, 3, 1, 1), pty = "s", mar = c(3, 4, 0, 0), tcl = 0.5, pin = c(3, 3)) X <- def$mono_rich Y <- def$resid P <- bquote(paste(P== .(sprintf("%1.2f", summary(lm.def1)[[1]][2,5])))) plot(jitter(X, amount = .1), Y, pch = 16, col = "gray", axes = T, xlim = c(-1, 30), ylim = YLIM, xlab = "Monoterpene richness", ylab = "Residuals from thickness-DBH relationship (mm)") # plot the points #axis(1) #axis(2, labels = NA) new <- data.frame(quantity = factor(2), mono_rich = seq(0, 30, by = 0.1)) pred <- predict(lm.def1, new, se.fit = T) lines(new$mono_rich, pred$fit) lines(new$mono_rich, pred$fit + pred$se.fit, lty = 2) lines(new$mono_rich, pred$fit - pred$se.fit, lty = 2) text(0, 17, P, pos = 4, cex = 1) mtext("Supplemental Figure 2", adj = 0, family = "Times", cex = 1.3) dev.off() ################################ ### SUPPLEMENTAL FIGURE 3 ### ################################ # PHYSIOLOGY # Bark thickness versus bark transpiration (Clement's data from Paracou) # try with species as covariate. First, eliminate any species with fewer than 3 indivs. lm.bark <- lm(Ecorce ~ DBH, dat = sap) summary(lm.bark) sap$resid[!is.na(sap$Ecorce) & !is.na(sap$DBH)] <- resid(lm.bark) lm1.i <- lm (RT.25 ~ resid, data = sap) lm2.i <- lme(RT.25 ~ resid, random = ~ 1 | name, data = sap, na.action = na.exclude) # controling for species diffs improves AIC #lm1.i <- lm (RT.25 ~ Ecorce, data = sap) #lm2.i <- lme(RT.25 ~ Ecorce, random = ~ 1 | name, data = sap, na.action = na.exclude) # controling for species diffs improves AIC #lm3.i <- lme(RT.25 ~ Ecorce * DBH, random = ~ 1 | name, data = sap, na.action = na.exclude) # adding DBH doesn't help at all! #lm3a.i <- lme(RT.25 ~ Ecorce, random = ~ DBH | name, data = sap, na.action = na.exclude, control = lmeControl(maxIter = 500, msMaxIter = 500)) # controling for Allowing random slopes with DBH also helps #lm3a.i <- lme(RT.25 ~ Ecorce, random = ~ 1 | DBH, data = sap, na.action = na.exclude, control = lmeControl(maxIter = 500, msMaxIter = 500)) # controling for Allowing random slopes with DBH also helps #lm3b.i <- lme(RT.25 ~ Ecorce, random = ~ DBH | Month/name, data = sap, na.action = na.exclude) # controling for Allowing random slopes with DBH also helps #lm3c.i <- lme(RT.25 ~ Ecorce, random = ~ DBH | Site/Month/name, data = sap, na.action = na.exclude) # controling for Allowing random slopes with DBH also helps #lm4.i <- lme(RT.25 ~ Ecorce, random = ~ 1 | Month, data = sap, na.action = na.exclude) # Block on Month of sampling #lm4a.i <- lme(RT.25 ~ Ecorce, random = ~ DBH | Month, data = sap, na.action = na.exclude, control = lmeControl(maxIter = 500, msMaxIter = 100)) # Won't converge AIC(lm1.i) AIC(lm2.i) AIC(lm3.i) AIC(lm3a.i) AIC(lm3b.i) AIC(lm3c.i) AIC(lm4.i) summary(lm2.i) # THis is the simplest and best-fitting model pdf(file = "figs/Supplemental Figure 3 2010 05 13.pdf", paper = "us", width = 8.5, height = 11) par(mfrow = c(1, 1), bty = "n", pin = c(3, 3), las = 1, tcl = 0.5) Xes <- data.frame(resid = seq(min(sap$resid, na.rm = T), max(sap$resid, na.rm = T), 0.1)) plot(sap$resid, sap$RT.25, ylab = expression(paste(CO[2], " efflux (", italic(ยต), "mol ", CO[2]%.%m^-2*s^-1, ")")), xlab = "Residuals from thickness-DBH relationsip (mm)", pch = 16, col = "gray", ylim = c(0, 5), xlim = c(-18, 18)) P <- bquote(paste(P== .(sprintf("%1.2f", summary(lm1.i)$coef[2, 4])))) pred <- predict(lm1.i, newdata = Xes, se.fit = T) lines(Xes$resid, pred$fit) lines(Xes$resid, pred$fit + pred$se.fit, lty = 2) lines(Xes$resid, pred$fit - pred$se.fit, lty = 2) text(-15, 4.25, P, pos = 4) mtext("Supplemental Figure 3", adj = 0, family = "Times", cex = 1.3, line = 2) dev.off()