In mdlama/milkweed:

library(milkweed)
library(dplyr)
library(ggplot2)
requirePackages(c("plot3D", "RColorBrewer", "gtable", "grid", "gridExtra"))
'%$%' <- magrittr::'%$%'
ipm <- mwIPM()

Let's prepare the data we will use for all vital rates.

x <- seq(from = 0, to = 160, length.out = 100)
y <- seq(from = 0, to = 6, length.out = 100)
M <- plot3D::mesh(x, y)

Flowering: Panel (a)

herb_avg <- ipm$data %>% 
  filter(!is.na(h_apical),
         !is.na(herb_avg),
         !is.na(fec.flower)) %$%
  herb_avg

Top

theme_set(theme_classic())

# Prepare data for heatmap
mydata <- tbl_df(data.frame(h_apical = as.vector(M$x),
                            herb_avg = as.vector(M$y),
                            prob.flower = predict(ipm$pars$flower.fit, 
                                                  newdata=data.frame(h_apical = as.vector(M$x),
                                                                     herb_avg = as.vector(M$y)),
                                                  type="Bertha")))

# Heatmap
p1 <- ggplot(mydata, aes(x = h_apical, 
                         y = herb_avg, 
                         z = prob.flower)) +
      geom_raster(aes(fill = prob.flower)) +
      geom_contour(colour = "white", alpha = 0.8) + 
      scale_fill_gradientn("Flowering\nProbability", 
                           colours=c("#00000000","#BBBBBBBB"),
                           limits=c(0, 1))

p1 <- p1 + theme(axis.line = element_blank()) + 
     scale_x_continuous(limits = c(0, max(M$x)), 
                        expand = c(0, 0)) + 
     scale_y_continuous(limits = c(0, 6), 
                        expand = c(0, 0)) +
     ylab("Herbivory severity")

# Add lines
herb_ex <- data.frame(yintercept = c(0, 
                                     mean(herb_avg), 
                                     max(herb_avg)))

p1 <- p1 + geom_hline(aes(yintercept = yintercept), 
                      data = herb_ex, 
                      linetype = c(1, 2, 5), 
                      col = RColorBrewer::brewer.pal(5, "YlOrRd")[2:4],
                      size=1)

# Move over a bit to match panel (b) below
p1 <- p1 + ggtitle("(a)") + 
           theme(axis.text.x = element_blank(),
                 axis.title.x = element_blank(),
                 axis.ticks.x = element_blank(),
                 plot.title = element_text(hjust = 0,
                                           margin = margin(b = 15, unit = "pt")))

p1

Bottom

min <- data.frame(h_apical = x, 
                 prob.flower = predict(ipm$pars$flower.fit, 
                                       newdata=data.frame(h_apical = x, 
                                                          herb_avg = herb_ex[1,1]), 
                                       type="Bertha"),
                 Herbivory = "min")
avg <- data.frame(h_apical = x, 
                  prob.flower = predict(ipm$pars$flower.fit, 
                                        newdata=data.frame(h_apical = x,
                                                           herb_avg = herb_ex[2,1]),
                                        type="Bertha"),
                  Herbivory = "avg")
max <- data.frame(h_apical = x, 
                  prob.flower = predict(ipm$pars$flower.fit, 
                                        newdata=data.frame(h_apical = x, 
                                                           herb_avg = herb_ex[3,1]),
                                        type="Bertha"),
                  Herbivory = "max")

mycurves <- rbind(min, avg, max) 
mycurves$Herbivory <- factor(mycurves$Herbivory, levels = c("min", "avg", "max"))

p2 <- ggplot(mycurves, aes(x = h_apical, 
                           y = prob.flower, 
                           col = Herbivory, 
                           linetype = Herbivory)) +
      geom_line(size = 1.0) + 
      scale_x_continuous(limits = c(0, max(x)), 
                         expand = c(0, 0)) + 
      scale_color_manual(values = RColorBrewer::brewer.pal(5, "YlOrRd")[2:4]) +
      # coord_cartesian(xlim = c(0, max(x))) + 
      xlab("Height (cm)") +
      ylab("Flowering Probability")     

p2 <- p2 + theme(legend.position = "none")

p2

Growth: Panel (b)

Load the aggregate data set & necessary packages

# Transform data
herb_avg <- ipm$data %>% 
  filter(!is.na(h_apical),
         !is.na(h_apical.next),
         !is.na(herb_avg),
         fec.flower == 1,
         surv == 1) %$%
  herb_avg

x <- seq(from = 0, to = 160, length.out = 100)
y <- seq(from = 0, to = 6, length.out = 100)
M <- plot3D::mesh(x, y)

Top

# Prepare data for heatmap
mydata <- tbl_df(data.frame(h_apical = as.vector(M$x),
                            herb_avg = as.vector(M$y),
                            h_apical.next = predict(ipm$pars$growth.fit, 
                                                    newdata=data.frame(h_apical = as.vector(M$x),
                                                                       herb_avg = as.vector(M$y)),
                                                    type="Bertha")))

# Heatmap
p3 <- ggplot(mydata, aes(x = h_apical, 
                         y = herb_avg, 
                         z = h_apical.next)) +
      geom_raster(aes(fill = h_apical.next)) +
      geom_contour(colour = "white", alpha = 0.8) + 
      scale_fill_gradientn("Height\n(Late Season)", 
                           colours=c("#00000000","#BBBBBBBB"),
                           limits=c(0, max(mydata$h_apical.next)))

# For some reason, coord_cartesian didn't work (left gap between tick marks and image)
p3 <- p3 + theme(axis.line = element_blank()) + 
     scale_x_continuous(limits = c(0, max(M$x)), 
                        expand = c(0, 0)) + 
     scale_y_continuous(limits = c(0, 6), 
                        expand = c(0, 0)) +
     ylab("Herbivory severity")

# Add lines
herb_ex <- data.frame(yintercept = c(0, 
                                     mean(herb_avg), 
                                     max(herb_avg)))

p3 <- p3 + geom_hline(aes(yintercept = yintercept), 
                      data = herb_ex, 
                      linetype = c(1, 2, 5), 
                      col = RColorBrewer::brewer.pal(5, "YlOrRd")[2:4],
                      size=1)

# Move over a bit to match panel (b) below
p3 <- p3 + ggtitle("(b)") + 
           theme(axis.text = element_blank(),
                 axis.title = element_blank(),
                 axis.ticks = element_blank(),
                 plot.title = element_text(hjust = 0,
                                           margin = margin(b = 15, unit = "pt")))

p3

Bottom

min <- data.frame(h_apical = x, 
                  h_apical.next = predict(ipm$pars$growth.fit, 
                                          newdata=data.frame(h_apical = x, 
                                                             herb_avg = herb_ex[1,1]), 
                                          type="Bertha"),
                  Herbivory = "min")
avg <- data.frame(h_apical = x, 
                  h_apical.next = predict(ipm$pars$growth.fit, 
                                          newdata=data.frame(h_apical = x,
                                                             herb_avg = herb_ex[2,1]),
                                          type="Bertha"),
                  Herbivory = "avg")
max <- data.frame(h_apical = x, 
                  h_apical.next = predict(ipm$pars$growth.fit, 
                                          newdata=data.frame(h_apical = x, 
                                                             herb_avg = herb_ex[3,1]),
                                          type="Bertha"),
                  Herbivory = "max")

mycurves <- rbind(min, avg, max) 
mycurves$Herbivory <- factor(mycurves$Herbivory, levels = c("min", "avg", "max"))

p4 <- ggplot(mycurves, aes(x = h_apical, 
                           y = h_apical.next, 
                           col = Herbivory, 
                           linetype = Herbivory)) +
      geom_line(size = 1.0) + 
      scale_x_continuous(limits = c(0, max(x)), 
                         expand = c(0, 0)) + 
      scale_color_manual(values = RColorBrewer::brewer.pal(5, "YlOrRd")[2:4]) +
      # coord_cartesian(xlim = c(0, max(x))) + 
      xlab("Height (cm)") +
      ylab("Height (Late Season)") +
      theme(legend.key.size =  unit(0.2, "in"),
            legend.position = c(1.0, 0.30))

p4

Now combine them!

g1 <- ggplotGrob(p1)
g2 <- ggplotGrob(p2)
g2 <- gtable::gtable_add_cols(g2, unit(0, "mm"))
g2 <- gtable::gtable_add_cols(g2, unit(0, "mm"))
g <- rbind(g1, g2, size="first")
g$widths <- grid::unit.pmax(g1$widths, g2$widths)

g3 <- ggplotGrob(p3)
g4 <- ggplotGrob(p4)
g4 <- gtable::gtable_add_cols(g4, unit(0, "mm"))
g4 <- gtable::gtable_add_cols(g4, unit(0, "mm"))
gg <- rbind(g3, g4, size="first")
gg$widths <- grid::unit.pmax(g3$widths, g4$widths)

p <- gridExtra::grid.arrange(g, gg, ncol=2)
# p <- gridExtra::grid.arrange(p2, p4, ncol=2)
ggsave("Figure2_FloweringAndGrowth.png", width=7, height=5, device = "png", p)