In mdlama/milkweed:

Initialization

library(milkweed)
library(dplyr)
library(ggplot2)
requirePackages(c("doParallel", "latex2exp", "gtable", "RColorBrewer", "gridExtra"))
library(doParallel) # Normally, I don't load suggested packages, but this one requires and loads parallel, foreach, and iterators,
                    #   so, you know, it's easier this way.
'%<>%' <- magrittr::'%<>%'
TeX <- latex2exp::TeX
ipm <- mwIPM()
mwCache <- getOption("milkweed.cache")

Analysis

Bootstrap elasticities

compute <- TRUE
parallel <- FALSE

if (!file.exists(file.path(mwCache,"figure6.RData")) | (compute)) {
  if (parallel) {
    numcl <- detectCores()-1
    registerDoParallel(cores = numcl)
    results <- tbl_df(foreach(i=1:20, .combine=bind_rows) %dopar% {
      library(milkweed)
      library(dplyr)
      '%<>%' <- magrittr::'%<>%'
      ipm %<>% bootIPM() %>% analyzeParameters(distpars = FALSE)
      thisone <- ipm$analysis$parameters %>% select(type, name, elasticity)
      thisone
    })
  } else {
    results <- tibble()
    for (i in 1:20) {
      ipm %<>% bootIPM() %>% analyzeParameters(distpars = FALSE)
      thisone <- ipm$analysis$parameters %>% select(type, name, elasticity)
      results <- bind_rows(results, thisone)
    }
  }

  save(results, file=file.path(mwCache,"figure6.RData"))
} else {
  load(file.path(mwCache,"figure6.RData"))
}

Spread elasticity over kernel for herbivory parameters

ipm %<>% analyzeStandard() %>% analyzeParameters(distpars = FALSE)
res <- ipm$analysis$parameters %>% 
    filter(type == "Pods") %>% 
    select(pars, elasticity)

pods_func <- function(x) {
  as.vector((ipm %>% setPodsMatrix(perturb = x))$kernels$K)
}

# Need 4 to handle all variables (even though only 2 with non-zero parameter values)
J <- numDeriv::jacobian(pods_func, rep(0,4))

# Loops!?  We don't need no stinkin loops...
S1 <- ipm$analysis$standard$sens*matrix(J[,1], nrow=ipm$N)*ipm$vars$h_apical$dx # Intercept
S2 <- ipm$analysis$standard$sens*matrix(J[,2], nrow=ipm$N)*ipm$vars$h_apical$dx # h_apical.next

E1 <- ipm$analysis$standard$elas*(matrix(J[,1], nrow=ipm$N)/(ipm$kernels$K/res$pars[1]))*ipm$vars$h_apical$dx
E2 <- ipm$analysis$standard$elas*(matrix(J[,2], nrow=ipm$N)/(ipm$kernels$K/res$pars[2]))*ipm$vars$h_apical$dx

# Sanity check - sum of Ei should equal calculated elasticities
abs(sum(sum(E1)) - res$elasticity[1])
abs(sum(sum(E2)) - res$elasticity[2])

Plotting

Caterpillar plot of elasticities

Prep the plot data

plotdata <- results %>% 
  tidyr::unite(label, type, name, sep="-", remove=FALSE) %>% 
  group_by(label) %>% 
  summarize(meanElasticity = mean(elasticity),
            n = n(),
            se = sd(elasticity)/sqrt(n),
            ci = qt(0.975, n-1)*se,
            type = last(type),
            name = last(name))

# Group parameters for plotting
new.levels <- c("Clonal", "Clonal", "Sexual", "Growth", "Herbivory", "Sexual", "Sexual", "Sexual", "Survival")
plotdata$newtype <- factor(new.levels[plotdata$type])

# Specify order of parameters in plot
newlabelorder <- c("Clonal-a", 
                   "Clonal-b", 
                   "Budlings-mean", 
                   "Budlings-sd", 
                   "Herbivory-pmunch", 
                   "Herbivory-mean", 
                   "Herbivory-sd", 
                   "Flowering-(Intercept)", 
                   "Flowering-h_apical", 
                   "Flowering-herb_avg", 
                   "Pods-(Intercept)", 
                   "Pods-h_apical.next", 
                   "Seedlings-mean", 
                   "Seedlings-sd", 
                   "Sexual-seedling.emergence", 
                   "Sexual-seeds.per.pod", 
                   "Growth-(Intercept)", 
                   "Growth-h_apical", 
                   "Growth-herb_avg", 
                   "Growth-sd", 
                   "Survival-(Intercept)",
                   "Survival-herb_avg")

# Corresponding mathematical symbols to above parameters
# These really ought to be done through plotmath notation:
#   https://stat.ethz.ch/R-manual/R-devel/library/grDevices/html/plotmath.html
fancypantslabels <- rev(newlabelorder)
fancypantslabels[1] <- TeX('$\\alpha$')
fancypantslabels[2] <- TeX('$\\beta_{z_{\\omega}}$')
fancypantslabels[3] <- TeX('$\\mu$')
fancypantslabels[4] <- TeX('$\\sigma$')
fancypantslabels[5] <- TeX('$p_{\\omega}$')
fancypantslabels[6] <- TeX('$\\mu$')
fancypantslabels[7] <- TeX('$\\sigma$')
fancypantslabels[8] <- TeX('$\\alpha$')
fancypantslabels[9] <- TeX('$\\beta_{z_{h}}$')
fancypantslabels[10] <- TeX('$\\beta_{z_{\\omega}}$')
fancypantslabels[11] <- TeX('$\\alpha$')
fancypantslabels[12] <- TeX('$\\beta_{z_{h}}$')
fancypantslabels[13] <- TeX('$\\beta_{z_{\\omega}}$')
fancypantslabels[14] <- TeX('$\\mu$')
fancypantslabels[15] <- TeX('$\\sigma$')
fancypantslabels[16] <- TeX('$\\alpha$')
fancypantslabels[17] <- TeX('$\\alpha$')
fancypantslabels[18] <- TeX('$\\alpha$')
fancypantslabels[19] <- TeX('$\\beta_{z_{h}}$')
fancypantslabels[20] <- TeX('$\\beta_{z_{\\omega}}$')
fancypantslabels[21] <- TeX('$\\sigma$')
fancypantslabels[22] <- TeX('$\\alpha$')
fancypantslabels[23] <- TeX('$\\beta_{z_{\\omega}}$')

# Legend ordering
newtypeorder <- c("Clonal", "Herbivory", "Sexual", "Growth", "Survival")

Render elasticity caterpillar plot

elastiplot <- plotdata %>% 
  mutate(label = factor(label, levels = rev(newlabelorder))) %>%
  mutate(newtype = factor(newtype, levels = newtypeorder)) %>%
  ggplot(aes(x = meanElasticity,
             y = label,
             color = newtype, 
             fill = newtype)) +
  geom_vline(xintercept = 0, 
             linetype = 1,
             size = 0.2) +
  geom_point() + 
  geom_errorbarh(aes(xmin = meanElasticity - ci,
                     xmax = meanElasticity + ci),
                height = 0.3)

# Highlight rects for functional forms
ym <- 0.25
al <- 0.1
r_survival <- annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1-ym, ymax = 2+ym, fill = "black", alpha = al)
r_growth <- annotate("rect", xmin = -Inf, xmax = Inf, ymin = 3-ym, ymax = 6+ym, fill = "black", alpha = al)
r_seedlings <- annotate("rect", xmin = -Inf, xmax = Inf, ymin = 9-ym, ymax = 10+ym, fill = "black", alpha = al)
r_pods <- annotate("rect", xmin = -Inf, xmax = Inf, ymin = 11-ym, ymax = 13+ym, fill = "black", alpha = al)
r_flowering <- annotate("rect", xmin = -Inf, xmax = Inf, ymin = 14-ym, ymax = 16+ym, fill = "black", alpha = al)
r_herbivory <- annotate("rect", xmin = -Inf, xmax = Inf, ymin = 17-ym, ymax = 19+ym, fill = "black", alpha = al)
r_budlings <- annotate("rect", xmin = -Inf, xmax = Inf, ymin = 20-ym, ymax = 21+ym, fill = "black", alpha = al)
r_bdlgs_per_stem <- annotate("rect", xmin = -Inf, xmax = Inf, ymin = 22-ym, ymax = 23+ym, fill = "black", alpha = al)

# Functional form labels
xc <- -0.4
sz <- 2.5
t_bdlgs_per_stem <- annotate("text", x = xc, y = 22.5, label = "Buds per stem", size = sz)
t_budlings <- annotate("text", x = xc, y = 20.5, label = "Budling recruits", size = sz)
t_herbivory <- annotate("text", x = xc, y = 18, label = "Herbivory", size = sz)
t_flowering <- annotate("text", x = xc, y = 15, label = "Flowering", size = sz)
t_pods <- annotate("text", x = xc, y = 12, label = "Pods", size = sz)
t_seedlings <- annotate("text", x = xc, y = 9.5, label = "Seedling recruits", size = sz)
t_seed_surv <- annotate("text", x = xc, y = 8, label = "Seed survival", size = sz)
t_seeds_per_pod <- annotate("text", x = xc, y = 7, label = "Seeds per pod", size = sz)
t_growth <- annotate("text", x = xc, y = 4.5, label = "Growth", size = sz)
t_survival <- annotate("text", x = xc, y = 1.5, label = "Ramet survival", size = sz)

# Throw those puppies in...
elastiplot <- elastiplot + 
  r_survival + 
  t_survival +
  r_growth + 
  t_growth +
  t_seeds_per_pod +
  t_seed_surv +
  r_seedlings + 
  t_seedlings +
  r_pods +
  t_pods +
  r_flowering +
  t_flowering +
  r_herbivory +
  t_herbivory +
  r_budlings +
  t_budlings +
  r_bdlgs_per_stem +
  t_bdlgs_per_stem

# Labels and such
elastiplot <- elastiplot + 
  scale_y_discrete(labels = rev(fancypantslabels)) +
  scale_x_continuous(limits = c(-1.2, 1.2)) +
  xlab("Parameter elasticity") +
  ylab("Parameters") + 
  theme_bw() +
  labs(color = "Parameter\nType",
       fill = "Parameter\nType") + 
  ggtitle("(a)") +
  theme(legend.background = element_rect(fill="lightgrey",
                                         size=0.1,
                                         linetype="solid"),
        legend.key.size = unit(0.18, "in"),
        legend.position = c(0.84, 0.25),
        plot.title = element_text(hjust = 0,
                                  margin = margin(b = 15, unit = "pt")))

elastiplot

Expanded herbivory elasticities

Prep data for heatmap of expanded elasticities

M <- plot3D::mesh(ipm$vars$h_apical$x, ipm$vars$h_apical$x)

plotdata <- tbl_df(data.frame(h_apical_t0 = as.vector(M$x),
                              h_apical_t1 = as.vector(M$y),
                              elasticity = as.vector(E1),
                              pname = "p[omega]"))
plotdata %<>% bind_rows(tbl_df(data.frame(h_apical_t0 = as.vector(M$x),
                                          h_apical_t1 = as.vector(M$y),
                                          elasticity = as.vector(E3),
                                          pname = "sigma")))

Render expanded elasticity plots for herbivory parameters

maxval <- max(abs(c(as.vector(E1), as.vector(E3))))
colorpal <- RColorBrewer::brewer.pal(3, "RdBu")

pmunch_elastiplot <- ggplot(plotdata, aes(x = h_apical_t0, 
                                          y = h_apical_t1, 
                                          z = elasticity)) +
  geom_raster(aes(fill = elasticity)) +
  facet_grid(pname ~ ., labeller=label_parsed) +
  geom_contour(colour = "black", alpha = 0.3) +
  scale_fill_gradient2(low = "blue",
                       mid = "white",
                       high = "red",
                       limits = c(-maxval,maxval))

pmunch_elastiplot <- pmunch_elastiplot + 
  theme(axis.line = element_blank()) + 
  scale_x_continuous(limits = c(min(ipm$vars$h_apical$x), 
                                max(ipm$vars$h_apical$x)), 
                     expand = c(0, 0)) + 
  scale_y_continuous(limits = c(min(ipm$vars$h_apical$x), 
                                max(ipm$vars$h_apical$x)), 
                     expand = c(0, 0)) + 
  xlab("Height at year t (cm)") + 
  ylab("Height at year t+1 (cm)") +
  labs(fill = "Elasticity") +
  theme_bw() +
  ggtitle("(b)") +
  theme(plot.title = element_text(hjust = 0,
                                  margin = margin(b = 15, unit = "pt")),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

pmunch_elastiplot

Put em together, it just makes sense!

g1 <- ggplotGrob(elastiplot)
g1 <- gtable::gtable_add_rows(g1, unit(0, "mm"))
g1 <- gtable::gtable_add_rows(g1, unit(0, "mm"))
g2 <- ggplotGrob(pmunch_elastiplot)
g <- cbind(g1, g2, size="first")

p <- gridExtra::grid.arrange(g1, g2, ncol=2)
ggsave("Figure6_ElasticityAnalysis.png", width=8, height=4, device = "png", p)