In mdlama/milkweed:

Initialization

library(milkweed)
library(dplyr)
library(ggplot2)
requirePackages(c("scales", "gridExtra"))
ipm <- mwIPM()

Plot Herbivory Distributions

Site

Make data for plot

N <- 1000
bx <- seq(from = 0, to = 6, length.out = N+1)
x <- 0.5*(bx[1:N] + bx[2:(N+1)])

site_ind <- c('BLD1','BLD2','GET','GRN','PWR','SKY','YTB','Bertha')
plotdata <- do.call(rbind, lapply(site_ind, 
                                  function (site) {
                                    data.frame(site = as.factor(site),
                                               x = x[2:(N-1)],
                                               y = ipm$pars$munched.fit[[site]]$predict(bx,
                                                                                        c(ipm$pars$munched.fit[[site]]$pmunch, 
                                                                                          ipm$pars$munched.fit[[site]]$fit$estimate),
                                                                                        justmunch=TRUE)[2:(N-1)])
                                  })
)

levels(plotdata$site)[levels(plotdata$site)=="Bertha"] <- "Combined"

Render the plot

plot_herb <- plotdata %>% ggplot(aes(x = x, y = y, fill = site, colour = site)) + 
  geom_line() + 
  geom_area(position = "identity", alpha = 0.3)

plot_herb <- plot_herb + theme_bw() +
  xlab("Herbivory severity") +
  ylab("Probability Density") +
  scale_x_log10(limits = c(-3, log10(6.1))) + 
  # scale_y_continuous(limits = c(0, 0.4)) +
  scale_fill_manual(values = c(scales::hue_pal()(7), NA)) + 
  scale_color_manual(values = c(scales::hue_pal()(7), "black")) +
  labs(colour="Sites", fill="Sites") +
  theme(legend.background = element_rect(fill="lightgrey",
                                         size=0.1,
                                         linetype="solid"),
        legend.key.size =  unit(0.18, "in"),
        legend.position = c(0.885, 0.72))

plot_herb <- plot_herb + ggtitle("(a)") + 
  theme(plot.title = element_text(hjust = 0, margin=margin(b = 15, unit = "pt")))

plot_herb

Year

Make data for plot

N <- 1000
bx <- seq(from = 0, to = 6, length.out = N+1)
x <- 0.5*(bx[1:N] + bx[2:(N+1)])

year_ind <- c('2013', '2014', '2015', '2016', '2017','Bertha')
plotdata <- do.call(rbind, lapply(year_ind, 
                                  function (year) {
                                    data.frame(year = as.factor(year),
                                               x = x[2:(N-1)],
                                               y = ipm$pars$munched.fit[[year]]$predict(bx,
                                                                                        c(ipm$pars$munched.fit[[year]]$pmunch, 
                                                                                          ipm$pars$munched.fit[[year]]$fit$estimate),
                                                                                        justmunch=TRUE)[2:(N-1)])
                                  })
)

levels(plotdata$year)[levels(plotdata$year)=="Bertha"] <- "Combined"

Render the plot

plot_herb <- plotdata %>% ggplot(aes(x = x, y = y, fill = year, colour = year)) + 
  geom_line() + 
  geom_area(position = "identity", alpha = 0.3)

plot_herb <- plot_herb + theme_bw() +
  xlab("Herbivory severity") +
  ylab("Probability Density") +
  scale_x_continuous(limits = c(0, 6)) + 
  # scale_y_continuous(limits = c(0, 0.4)) +
  scale_fill_manual(values = c(scales::hue_pal()(5), NA)) + 
  scale_color_manual(values = c(scales::hue_pal()(5), "black")) +
  labs(colour="Years", fill="Years") +
  theme(legend.background = element_rect(fill="lightgrey",
                                         size=0.1,
                                         linetype="solid"),
        legend.key.size =  unit(0.18, "in"),
        legend.position = c(0.885, 0.72))

plot_herb <- plot_herb + ggtitle("(a)") + 
  theme(plot.title = element_text(hjust = 0, margin=margin(b = 15, unit = "pt")))

plot_herb

Plot Budling and Seedling Distributions

Make data for plot

N <- 1000
bx <- seq(from = 0, to = 160, length.out = N+1)
x <- 0.5*(bx[1:N] + bx[2:(N+1)])

# First add budling distributions
site_ind <- c('BLD1','BLD2','GET','GRN','PWR','SKY','YTB','Bertha')
plotdata <- do.call(rbind, lapply(site_ind, 
                                  function (site) {
                                    data.frame(site = site,
                                               x = x[3:N-1],
                                               y = ipm$pars$budling.fit[[site]]$predict(bx,
                                                                                        ipm$pars$budling.fit[[site]]$fit$estimate)[3:N-1])
                                  })
)

# Now add seedling distributions
plotdata <- bind_rows(plotdata,
                      data.frame(site = "Seedlings",
                                 x = x[3:N-1],
                                 y = ipm$pars$seedling.fit$predict(bx, 
                                                                   ipm$pars$seedling.fit$fit$estimate)[3:N-1]))

plotdata$site <- factor(plotdata$site, levels = c('BLD1','BLD2','GET','GRN','PWR','SKY','YTB','Bertha','Seedlings'))

levels(plotdata$site)[levels(plotdata$site)=="Bertha"] <- "Combined"

#plotdata$site <- factor(plotdata$site, levels = c('BLD1','BLD2','PWR','SKY','YTB','Combined','Seedlings'))

Render the plot

plot_recr <- plotdata %>% ggplot(aes(x = x, 
                                     y = y, 
                                     fill = site, 
                                     colour = site,
                                     linetype = site)) + 
  geom_line() + 
  geom_area(position = "identity", alpha = 0.3)

plot_recr <- plot_recr + theme_bw() +
  xlab("Height (cm)") +
  ylab("Probability Density") +
  scale_x_continuous(limits = c(0, 160)) +
  scale_fill_manual(values = c(scales::hue_pal()(7), NA, NA)) + 
  scale_color_manual(values = c(scales::hue_pal()(7), "black", "black")) + 
  scale_linetype_manual(values = c(1,1,1,1,1,1,1,1,2)) +
  labs(colour="Sites", fill="Sites", linetype="Sites") +
  theme(legend.background = element_rect(fill="lightgrey",
                                         size=0.1,
                                         linetype="solid"),
        legend.key.size =  unit(0.2, "in"),
        legend.position = c(0.875, 0.65))

plot_recr <- plot_recr + ggtitle("(b)") + 
  theme(plot.title = element_text(hjust = 0, margin=margin(b = 15, unit = "pt")))

plot_recr

Render and save final plot

pall <- gridExtra::grid.arrange(plot_herb, plot_recr, ncol=2)
ggsave("AppendixFigure2_RecruitAndHerbDist.png", height=4, width=12, device = "png", pall)