#=============================================================================================================# #### R code for "Prairie fires improve the mating scence and reproduction of Echinacea angustifolia" #### #=============================================================================================================# #### 1. Functions #### require(ggplot2) require(grid) require(gridExtra) #=============================================================================================================# #### 2. Read in Echinacea data #### data <- read.csv("http://echinaceaproject.org/data/echinacea_sppFireData.csv") # data <- read.csv("forWebsite/echinacea_sppFireData.csv") str(data) #=============================================================================================================# #### 3. Metadata #### # $ unit: burn unit at Staffanson Prairie Preserve, "east" or "west" # $ year: year in which reproductive data, 1996-2014 # $ burned: indicates whether unit was burned prior to the growing season, "burned" or "unburned" # $ flaPlas: count of flowering plants in burn unit # $ totalPlas: count of total living plants, flowering plants + nonflowering plants # $ flw: proportion of plants in flower # $ hdCt: average number of flower heads per flowering plant # $ flwDensity: flowering plants per hectare, count of flowering plants divided by unit area # $ nn1.hmean: harmonic mean of distance to nearest flowering neighbor, in meters # $ nn2.hmean: harmonic mean of distance to second nearest flowering neighbor, in # $ nn3.hmean: harmonic mean of distance to third nearest flowering neighbor, in meters # $ nn4.hmean: harmonic mean of distance to fourth nearest flowering neighbor, in meters # $ meanSdSet: average proportion of achenes which contain fertilized embryos # $ meanAcheneCt: average number of achenes produced per flower head, includes fertilized and unfertilized achenes # $ meanFecundity: average number of fertilized achenes produced per plant, product of individual seed set, achene count and head count # $ meanSp1: average number of days styles persist on a flower head # $ mean.sync: average pairwise temporal flowering synchrony, proportion of overlapping flowering days # $ sdsd: standard deviation of flowering start date # $ mean.duration: average flowering duration, in days # $ lwr95.sp1: lower bound of bootstrapped 95% confidence interval for mean style persistence # $ up95.sp1: upper bound of bootstrapped 95% confidence interval for mean style persistence # $ lwr95.sdSet: lower bound of bootstrapped 95% confidence interval for mean seed set # $ up95.sdSet: upper bound of bootstrapped 95% confidence interval for mean seed set #=============================================================================================================# #### 4. Analyze relationship between mating potential and reproductive success #### #### > Synchrony vs. seed set #### plot(meanSdSet~mean.sync, data) lm.sync.ss <- lm(meanSdSet~mean.sync, data) lm.sync.ss.remove1997 = lm(meanSdSet~mean.sync, data[-2,]) abline(lm.sync.ss) abline(lm.sync.ss.remove1997, col = "red") # plot(lm.sync.ss, which = 1) # plot(lm.sync.ss, which = 2) # plot(lm.sync.ss.remove1997, which = 1) # plot(lm.sync.ss.remove1997, which = 2) summary(lm.sync.ss) summary(lm.sync.ss.remove1997) #### > Isolation vs. seed set #### ### flowering plant density plot(meanSdSet~log(flwDensity), data) lm.dens.ss <- lm(meanSdSet~log(flwDensity), data) # plot(lm.dens.ss, which = 1) # plot(lm.dens.ss, which = 2) summary(lm.dens.ss) #### > harmonic mean of distance to nearest flowering neighbor #### plot(meanSdSet~log(nn1.hmean), data) lm.nn1.ss <- lm(meanSdSet~log(nn1.hmean), data) # plot(lm.nn1.ss, which = 1) # plot(lm.nn1.ss, which = 2) summary(lm.nn1.ss) #### > harmonic mean of distance to second nearest flowering neighbor #### plot(meanSdSet~log(nn2.hmean), data) lm.nn2.ss <- lm(meanSdSet~log(nn2.hmean), data) # plot(lm.nn2.ss, which = 1) # plot(lm.nn2.ss, which = 2) summary(lm.nn2.ss) #### > harmonic mean of distance to third nearest flowering neighbor #### plot(meanSdSet~log(nn3.hmean), data) lm.nn3.ss <- lm(meanSdSet~log(nn3.hmean), data) # plot(lm.nn3.ss, which = 1) # plot(lm.nn3.ss, which = 2) summary(lm.nn3.ss) #### > harmonic mean of distance to fourth nearest flowering neighbor #### plot(meanSdSet~log(nn4.hmean), data) lm.nn4.ss <- lm(meanSdSet~log(nn4.hmean), data) # plot(lm.nn4.ss, which = 1) # plot(lm.nn4.ss, which = 2) summary(lm.nn4.ss) #### > seed set vs. style persistence #### plot(meanSdSet~meanSp1, data) lm.ss.sp <- lm(meanSdSet~meanSp1, data) plot(lm.ss.sp, which = 1) plot(lm.ss.sp, which = 2) summary(lm.ss.sp) #=============================================================================================================# #### 5. Compare Echinacea reproduction in burned and unburned units #### ### use Wilcoxon rank-sum tests to test for differences between burned and unburned units #### > Annual reproduction #### boxplot(flw~burned, data, ylab = "Proportion flowering") wilcox.test(flw~burned, data) boxplot(hdCt~burned, data, ylab = "Mean head count") wilcox.test(hdCt~burned, data) boxplot(meanSdSet~burned, data, ylab = "Mean seed set") wilcox.test(meanSdSet~burned, data) boxplot(meanAcheneCt~burned, data, ylab = "Mean achene count") wilcox.test(meanAcheneCt~burned, data) boxplot(meanFecundity~burned, data, ylab = "Mean fecundity") wilcox.test(meanFecundity~burned, data) boxplot(meanSp1~burned, data, ylab = "Mean style persistence (days)") wilcox.test(meanSp1~burned, data) #### > 5b. Spatial measures of mating poential #### boxplot(flwDensity~burned, data, ylab = "Mean flowering density (plants/ha)") wilcox.test(flwDensity~burned, data) boxplot(nn1.hmean~burned, data, ylab = "Harmonic mean of distance to first nearest neighbor (m)") wilcox.test(nn1.hmean~burned, data) boxplot(nn2.hmean~burned, data, ylab = "Harmonic mean of distance to second nearest neighbor (m)") wilcox.test(nn2.hmean~burned, data) boxplot(nn3.hmean~burned, data, ylab = "Harmonic mean of distance to third nearest neighbor (m)") wilcox.test(nn3.hmean~burned, data) boxplot(nn4.hmean~burned, data, ylab = "Harmonic mean of distance to fourth nearest neighbor (m)") wilcox.test(nn4.hmean~burned, data) #### > 5c. Temporal measures of mating potential #### boxplot(mean.sync ~ burned, data) wilcox.test(mean.sync ~ burned, data) boxplot(sdsd~burned, data, ylab = "Standard deviation of flowering start day") wilcox.test(sdsd~burned, data) boxplot(mean.duration~burned, data, ylab = "Mean flowering duration (days)") wilcox.test(mean.duration~burned, data) #=============================================================================================================# #=============================================================================================================# #### 6. Figure 2 #### predict = data.frame(nn2.hmean = seq(1.5,16, 0.5)) predict$mean.sync = seq(0.3072, 0.8036, 0.4964/29) predict$meanSp1 = seq(1.317, 1.714, 0.397/29) predict$lm.sync.ss.predict = predict(lm.sync.ss, predict, type = "response") predict$lm.sync.ss.predict.remove1997 = predict(lm.sync.ss.remove1997, predict, type = "response") predict$lm.nn2.ss.predict = predict(lm.nn2.ss, predict, type = "response") predict$lm.sp.ss.predict = predict(lm.ss.sp, predict, type = "response") ptsize = 1.5 pt2size = 5 letsize = 4 ciColor = "black" baseSize = 9 fig2a <- ggplot(data[!data$meanSdSet %in% NA, ], aes(x=mean.sync, y = meanSdSet, shape = unit, fill = unit)) + geom_point(data = data[data$burned %in% "burned" & !data$meanSdSet %in% NA,], aes(x = mean.sync, y = meanSdSet), colour="red", size=pt2size, alpha=.75, inherit.aes=FALSE) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + ## Asterisk for 1997 geom_point(data = data[data$mean.sync > 0.8,], aes(x=mean.sync +0.025, y = meanSdSet+0.025), shape = 8, size=ptsize-0.5, color = "grey20", inherit.aes = F) + geom_line(data = predict, aes(x=mean.sync, y = lm.sync.ss.predict), linetype = 2, inherit.aes=FALSE) + geom_line(data = predict, aes(x=mean.sync, y = lm.sync.ss.predict.remove1997), linetype = 1, inherit.aes=FALSE) + annotate("text", 0.27, 0.8, size = letsize, hjust = 0, vjust = 1, label = "A") + annotate("text", 0.36, 0, size = 3, hjust = 0, vjust = 0, label = "N = 20", fontface = 3) + scale_shape_manual(values=c(24,21)) + scale_fill_manual(values=c("black","white")) + scale_x_continuous("Mean pairwise\nflowering synchrony", limits=c(0.27,0.9), breaks = seq(0.3,0.75,0.15)) + scale_y_continuous("Mean seed set", limits=c(0,0.8), breaks=seq(0,0.8,.2), labels = c("0.00", "0.20", "0.40", "0.60", "0.80")) + theme_classic(base_size = baseSize) + theme(legend.position = "none", axis.text.y = element_text(hjust = 0.5), axis.title.y = element_text(vjust = 1, hjust = .5, size = 10), plot.margin = margin(r = -0.3, unit = "cm")) fig2b <- ggplot(data[!data$meanSdSet %in% NA,], aes(x=nn2.hmean, y = meanSdSet, shape = unit, fill = unit)) + geom_point(data = data[data$burned %in% "burned" & !data$meanSdSet %in% NA,], aes(x = nn2.hmean, y = meanSdSet), colour="red", size=pt2size, alpha=.75, inherit.aes=FALSE) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + ## Asterisk for 1997 geom_point(data = data[data$mean.sync > 0.8,], aes(x=nn2.hmean +1, y = meanSdSet+0.025), shape = 8, color = "grey20", size=ptsize-0.5, inherit.aes = F) + geom_line(data = predict, aes(x=nn2.hmean, y = lm.nn2.ss.predict), linetype = 1, inherit.aes=FALSE) + annotate("text", 1.5, 0.8, size = letsize, hjust = 0, vjust = 1, label = "B") + annotate("text", 2, 0, size = 3, hjust = 0, vjust = 0, label = "N = 20", fontface = 3) + scale_shape_manual(values=c(24,21), name = "Burn Unit", labels = c("East", "West")) + scale_fill_manual(values=c("black","white"), name = "Burn Unit", labels = c("East", "West")) + scale_x_log10("Mean distance to second\nnearest flowering plant (m)", limits = c(1.5,20), breaks = c(2,4,8,16)) + scale_y_continuous("Mean seed set", limits=c(0,0.8), breaks=seq(from = 0, to = 0.8, by =.2)) + theme_classic(base_size = baseSize) + theme(legend.position = "none", axis.text.y = element_text(colour = "white"), axis.title.y = element_text(vjust = 1, hjust = .5, colour = "white"), plot.margin = margin(l = -0.15, r=-0.15, unit = "cm")) fig2c <- ggplot(data[!data$meanSp1 %in% NA,], aes(x=meanSp1, y = meanSdSet, shape = unit, fill = unit)) + geom_point(data = data[data$burned %in% "burned" & !data$meanSdSet %in% NA,], aes(x = meanSp1, y = meanSdSet), colour="red", size=pt2size, alpha=.75, inherit.aes=FALSE) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + ## Asterisk for 1997 geom_point(data = data[data$mean.sync > 0.8,], aes(x=meanSp1 +0.1, y = meanSdSet+0.025), shape = 8, size=ptsize-0.5, color = "grey20", inherit.aes =F) + stat_smooth(data = data, aes(x=meanSp1, y = meanSdSet), method = "lm", se = F, colour = "black", size = 0.5, inherit.aes = F) + annotate("text", 1.25, 0.8, size = letsize, hjust = 0, vjust = 1, label = "C") + annotate("text", 1.5, 0, size = 3, hjust = 0, vjust = 0, label = "N = 16", fontface = 3) + scale_shape_manual(values=c(24,21), name = "Burn Unit", labels = c("East", "West")) + scale_fill_manual(values=c("black","white"), name = "Burn Unit", labels = c("East", "West")) + scale_x_continuous("Mean style\npersistence (days)", limits = c(1.25,3.25), breaks = seq(1.5,3,0.5)) + scale_y_continuous("Mean seed set", limits=c(0,0.8), breaks=seq(from = 0, to = 0.8, by =.2)) + geom_point(data = NULL, aes(x = 2.65, y = 0.725), colour="red", size=5, alpha=.75, inherit.aes=FALSE) + geom_point(data = NULL, aes(x = 3, y = 0.725), colour="red", size=5, alpha=.75, inherit.aes=FALSE) + geom_point(data = NULL, aes(x=3, y = 0.725), shape = 21, fill = "white", size=1.5) + geom_point(data = NULL, aes(x=3, y = 0.65), shape = 21, fill = "white", size=1.5) + geom_point(data = NULL, aes(x=2.65, y = 0.65), shape = 24, fill = "black", size=1.5) + geom_point(data = NULL, aes(x=2.65, y = 0.725), shape = 24, fill = "black", size=1.5) + annotate("text", 2.65, 0.8, size = 2.5, hjust = 0.5, vjust = 0.5, label = "East") + annotate("text", 3, 0.8, size = 2.5, hjust = 0.5, vjust = 0.5, label = "West") + annotate("text", 2.4, 0.725, size = 2.5, hjust = 1, vjust = 0.5, label = "Burned") + annotate("text", 2.4, 0.65, size = 2.5, hjust = 1, vjust = 0.5, label = "Unburned") + theme_classic(base_size = baseSize) + theme(legend.position = "none", axis.text.y = element_text(colour = "white"), axis.title.y = element_text(vjust = 1, hjust = .5, colour = "white"), plot.margin = margin(l = -0.3, unit = "cm")) grid.arrange(fig2a, fig2b, fig2c, ncol = 3) #=============================================================================================================# #### 7. Figure 3 #### ptsize = 2 pt2size = 7 letsize = 4 baseText = 8 axisTitle = 8 xlabels = c(1996, "", 1998, "", 2000, "", 2002, "", 2004, "", 2006, "", 2008, "", 2010, "", 2012, "", 2014) data$year2 <- as.factor(data$year) fig3a <- ggplot(data=data[!data$mean.sync %in% NA,], aes(x=year2, y=mean.sync, group=unit)) + geom_point(data = data[data$burned %in% "burned" & !data$mean.sync %in% NA,], aes(x = year2, y = mean.sync), colour="red", size=pt2size, inherit.aes=FALSE) + geom_line(data = data[!data$mean.sync %in% NA & data$unit %in% "east",], aes(x=year2, y=mean.sync, group = unit), linetype = 1, inherit.aes=F) + geom_line(data = data[!data$mean.sync %in% NA & data$unit %in% "west",], aes(x=year2, y=mean.sync, group = unit), linetype = 2, inherit.aes=F) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + ## Asterisk for 1997 geom_point(data = data[data$mean.sync > 0.8,], aes(x=year2, y = mean.sync+0.09), shape = 8, color = "grey20", size=ptsize-0.5, inherit.aes = F) + annotate("text", 0.75, 0.9, size = letsize, hjust = 0, vjust = 1, label = "A") + scale_x_discrete("") + scale_y_continuous("Mean pairwise\nflowering synchrony", limits=c(0.25,0.9), breaks = seq(0.3,0.9,0.15)) + scale_shape_manual(values=c(24,21), name = "Burn Unit", labels = c("East", "West")) + scale_fill_manual(values=c("black","white"), name = "Burn Unit", labels = c("East", "West")) + scale_linetype_manual(values=c("solid","dashed")) + theme_classic(base_size = baseText) + theme(legend.position = "none", plot.margin=unit(c(0.5,0.5,-0.125, 0.5), "cm"), legend.justification = c(1,0.1), legend.direction = "horizontal", axis.text.x = element_blank(), axis.title.x = element_text(vjust = .5, hjust = .5, size = axisTitle), axis.title.y = element_text(vjust = 1, hjust = .5, size = axisTitle)) fig3b <- ggplot(data=data[!data$sdsd %in% NA, ], aes(x=year2, y=sdsd, group=unit)) + geom_point(data = data[data$burned %in% "burned" & !data$sdsd %in% NA,], aes(x = year2, y = sdsd), colour="red", size=pt2size, inherit.aes=FALSE) + geom_line(aes(linetype=unit)) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + ## Asterisk for 1997 geom_point(data = data[data$mean.sync > 0.8,], aes(x=year2, y = sdsd+0.85), shape = 8, color = "grey20", size=ptsize-0.5, inherit.aes = F) + annotate("text", 0.75, 10, size = letsize, hjust = 0, vjust = 1, label = "B") + scale_x_discrete("") + scale_y_continuous("Standard deviation\nof start day", limits=c(2,10), breaks=seq(2,10,by=2), labels = c("0.0", "1.0", "2.0", "3.0", "4.0")) + scale_shape_manual(values=c(24,21)) + scale_fill_manual(values=c("black","white")) + scale_linetype_manual(values=c("solid","dashed")) + theme_classic(base_size = baseText) + theme(legend.position = "none", plot.margin=unit(c(0,0.5,0,0.5), "cm"), axis.text.x = element_blank(), axis.title.x = element_text(vjust = .5, hjust = .5, size = axisTitle), axis.title.y = element_text(vjust = 1, hjust = .5, size = axisTitle)) fig3c <- ggplot(data=data[!data$mean.duration %in% NA, ], aes(x=year2, y=mean.duration, group=unit)) + geom_point(data = data[data$burned %in% "burned" & !data$mean.duration %in% NA,], aes(x = year2, y = mean.duration), colour="red", size=pt2size, inherit.aes=FALSE) + geom_line(aes(linetype=unit)) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + ## Asterisk for 1997 geom_point(data = data[data$mean.sync > 0.8,], aes(x=year2, y = mean.duration+2), shape = 8, color = "grey20", size=ptsize-0.5, inherit.aes = F) + annotate("text", 0.75, 24, size = letsize, hjust = 0, vjust = 1, label = "C") + scale_x_discrete("Year") + scale_y_continuous("Mean flowering\nduration (days)", limits=c(11,25), breaks = c(12, 16, 20,24), labels=c(" 12", " 16", " 20", " 24")) + scale_shape_manual(values=c(24,21)) + scale_fill_manual(values=c("black","white")) + scale_linetype_manual(values=c("solid","dashed")) + geom_point(data = NULL, aes(x = 9.125, y = 23), colour="red", size=5, alpha=.75, inherit.aes=FALSE) + geom_point(data = NULL, aes(x = 10.125, y = 23), colour="red", size=5, alpha=.75, inherit.aes=FALSE) + geom_point(data = NULL, aes(x=9.125, y = 20.5), shape = 24, fill = "black", size=1.5) + geom_point(data = NULL, aes(x=9.125, y = 23), shape = 24, fill = "black", size=1.5) + geom_point(data = NULL, aes(x=10.125, y = 20.5), shape = 21, fill = "white", size=1.5) + geom_point(data = NULL, aes(x=10.125, y = 23), shape = 21, fill = "white", size=1.5) + annotate("text", 9.125, 25, size = 2.5, hjust = 0.5, vjust = 0.5, label = "East") + annotate("text", 10.125, 25, size = 2.5, hjust = 0.5, vjust = 0.5, label = "West") + annotate("text", 8.625, 23, size = 2.5, hjust = 1, vjust = 0.5, label = "Burned") + annotate("text", 8.625, 20.5, size = 2.5, hjust = 1, vjust = 0.5, label = "Unburned") + theme_classic(base_size = baseText) + theme(legend.position = "none", plot.margin=unit(c(-0.125,0.5,0.25,0.5), "cm"), axis.title.x = element_text(vjust = .5, hjust = .5, size = axisTitle), axis.title.y = element_text(vjust = 1, hjust = .5, size = axisTitle)) grid.arrange(fig3a, fig3b, fig3c, ncol = 1) #=============================================================================================================# #### 8. Figure 4 #### baseText = 8 ptsize = 1.5 pt2size = 5 letsize = 4 ciColor = "black" baseSize = 9 axisTitle = 8 xlabels = c(1996, "", 1998, "", 2000, "", 2002, "", 2004, "", 2006, "", 2008, "", 2010, "", 2012, "", 2014, "", 2016) fig4a <- ggplot(data=data, aes(x=year, y=nn2.hmean, group=unit)) + geom_point(data = data[data$burned %in% "burned",], aes(x = year, y = nn2.hmean), colour="red", size=pt2size, inherit.aes=FALSE) + geom_line(data = data[data$unit %in% "east",], aes(x=year, y=nn2.hmean, group = unit), linetype = 1, inherit.aes=F) + geom_line(data = data[data$unit %in% "west",], aes(x=year, y=nn2.hmean, group = unit), linetype = 2, inherit.aes=F) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + annotate("text", 1995.875, 32, size = letsize, hjust = 1, vjust = 1, label = "A") + scale_x_continuous("", limits=c(1995.5,2016), breaks=seq(1996,2016, by=1), labels = xlabels) + scale_y_continuous("Mean distance to second \nnearest neighbor (m)", trans = "log", limits=c(1,32), breaks=c(1,2,4,8,16, 32), labels = c(" 1", " 2", " 4", " 8", " 16", "32")) + scale_shape_manual(values=c(24,21), name = "Burn Unit", labels = c("East", "West")) + scale_fill_manual(values=c("black","white"), name = "Burn Unit", labels = c("East", "West")) + scale_linetype_manual(values=c("solid","dashed"), name = "Burn Unit", labels = c("East", "West")) + geom_point(data = NULL, aes(x = 2013, y = 20), colour="red", size=5, alpha=.75, inherit.aes=FALSE) + geom_point(data = NULL, aes(x = 2015.5, y = 20), colour="red", size=5, alpha=.75, inherit.aes=FALSE) + geom_point(data = NULL, aes(x=2013, y = 20), shape = 24, fill = "black", size=1.5) + geom_point(data = NULL, aes(x=2013, y = 11), shape = 24, fill = "black", size=1.5) + geom_point(data = NULL, aes(x=2015.5, y = 11), shape = 21, fill = "white", size=1.5) + geom_point(data = NULL, aes(x=2015.5, y = 20), shape = 21, fill = "white", size=1.5) + annotate("text", 2013, 32, size = 2.5, hjust = 0.5, vjust = 0.5, label = "East") + annotate("text", 2015.5, 32, size = 2.5, hjust = 0.5, vjust = 0.5, label = "West") + annotate("text", 2012, 20, size = 2.5, hjust = 1, vjust = 0.5, label = "Burned") + annotate("text", 2012, 11, size = 2.5, hjust = 1, vjust = 0.5, label = "Unburned") + theme_classic(base_size = baseText) + theme(legend.position = "none", plot.margin=unit(c(0.125,0.5,-0.25, 0.5), "cm"), legend.direction = "horizontal", axis.text.x = element_blank(), axis.title.x = element_text(vjust = .5, hjust = .5, size = axisTitle), axis.title.y = element_text(vjust = 1, hjust = .5, size = axisTitle)) fig4b <- ggplot(data=data, aes(x=year, y=flw, group=unit)) + geom_point(data = data[data$burned %in% "burned",], aes(x = year, y = flw), colour="red", size=pt2size, inherit.aes=FALSE) + geom_line(aes(linetype=unit)) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + annotate("text", 1995.875, 0.82, size = letsize, hjust = 1, vjust = 1, label = "B") + scale_x_continuous("Year", limits=c(1995.5,2016), breaks=seq(1996,2016, by=1), labels = xlabels) + scale_y_continuous("Proportion of living\nplants in flower", limits=c(0,0.82), breaks=seq(0,0.8,by=.2)) + scale_shape_manual(values=c(24,21)) + scale_fill_manual(values=c("black","white")) + scale_linetype_manual(values=c("solid","dashed")) + theme_classic(base_size = baseText) + theme(legend.position = "none", plot.margin=unit(c(0.0,0.5,0.3125,0.5), "cm"), axis.title.x = element_blank(), axis.text.x = element_blank(), axis.title.y = element_text(vjust = 1, hjust = .5, size = axisTitle)) fig4c <- ggplot(data=data, aes(x=year, y=hdCt, group=unit)) + geom_point(data = data[data$burned %in% "burned",], aes(x = year, y = hdCt), colour="red", size=pt2size, inherit.aes=FALSE) + geom_line(aes(linetype=unit)) + geom_point(aes(shape=unit, fill=unit), size=ptsize) + annotate("text", 1995.875, 3.5, size = letsize, hjust = 1, vjust = 1, label = "C") + scale_x_continuous("Year", limits=c(1995.5,2016), breaks=seq(1996,2016, by=1), labels = xlabels) + scale_y_continuous("Mean head count\nper flowering plant", limits=c(0.5,3.5), breaks=seq(0,3.5,by=0.5)) + scale_shape_manual(values=c(24,21)) + scale_fill_manual(values=c("black","white")) + scale_linetype_manual(values=c("solid","dashed")) + theme_classic(base_size = baseText) + theme(legend.position = "none", plot.margin=unit(c(-0.125,0.5,0.125,0.5), "cm"), axis.title.x = element_text(vjust = .5, hjust = .5, size = axisTitle), axis.title.y = element_text(vjust = 1, hjust = .5, size = axisTitle)) grid.arrange(fig4a, fig4b, fig4c, ncol = 1)