In mdlama/milkweed:

Initialization

source("../../../.Rprofile")

Linear

ipm <- mwIPM(list(compute = TRUE,
                  mdlargs = list(method = "linear",
                                 input = "full")))

sites <- c("Bertha", "BLD1", "BLD2", "PWR", "SKY", "YTB")
sitedata <- bind_cols(tbl_df(t(sapply(ipm$pars$munched.fit[-1], function (x) {c(x$fit$estimate)}))),
                      tbl_df(sapply(ipm$pars$munched.fit[-1], function (x) {x$pmunch})),
                      tbl_df(data.frame(site = sites[-1])))

meanlog <- seq(from = -2.25, 
               to = 0.0, 
               length.out = 50)
sdlog <- seq(from = 0.9*min(sitedata$sdlog), 
             to = 1.1*max(sitedata$sdlog), 
             length.out = 50)
M <- mesh(meanlog, sdlog)

Compute pop growth as function of mean and sd (linear)

compute <- FALSE

if (!file.exists(mwROOT("data","calculated","figure7a.RData")) | (compute)) {
  plotdata <- tbl_df(
    data.frame(
      meanlog = as.vector(M$x),
      sdlog = as.vector(M$y),
      type = "linear"
    )
  )

  pmunch <- ipm$pars$munched.fit[["Bertha"]]$pmunch
  munched <- ipm$pars$munched.fit[["Bertha"]]$predict
  b <- ipm$vars$log_herb_avg$b
  plotdata$growth_rate = sapply(1:nrow(plotdata), 
                                function (x) { 
                                  ipm %>% 
                                    setHerbivoryMatrix(dist.herb = munched(b, 
                                                                           pars = c(pmunch, 
                                                                                    plotdata$meanlog[x],
                                                                                    plotdata$sdlog[x]))) %>% 
                                    analyzeGrowthRate()
                                })

  save(plotdata, file=mwROOT("data","calculated","figure7a.RData"))
} else {
  load(mwROOT("data","calculated","figure7a.RData"))
}

Exponential

ipm <- mwIPM(list(compute = TRUE,
                  mdlargs = list(method = "exp",
                                 input = "full")))

sites <- c("Bertha", "BLD1", "BLD2", "PWR", "SKY", "YTB")
sitenames <- c("Combined", "BLD1", "BLD2", "PWR", "SKY", "YTB")
sitedata <- bind_cols(tbl_df(t(sapply(ipm$pars$munched.fit, function (x) {c(x$fit$estimate)}))),
                      tbl_df(sapply(ipm$pars$munched.fit, function (x) {x$pmunch})),
                      tbl_df(data.frame(site = sitenames)))

Compute pop growth as function of mean and sd (exponential)

compute <- FALSE

if (!file.exists(mwROOT("data","calculated","figure7b.RData")) | (compute)) {
  plotdatab <- tbl_df(
    data.frame(
      meanlog = as.vector(M$x),
      sdlog = as.vector(M$y),
      type = "exp"
    )
  )

  pmunch <- ipm$pars$munched.fit[["Bertha"]]$pmunch
  munched <- ipm$pars$munched.fit[["Bertha"]]$predict
  b <- ipm$vars$log_herb_avg$b
  plotdatab$growth_rate = sapply(1:nrow(plotdatab), 
                                function (x) { 
                                  ipm %>% 
                                    setHerbivoryMatrix(dist.herb = munched(b, 
                                                                           pars = c(pmunch, 
                                                                                    plotdata$meanlog[x],
                                                                                    plotdata$sdlog[x]))) %>% 
                                    analyzeGrowthRate()
                                })

  save(plotdatab, file=mwROOT("data","calculated","figure7b.RData"))
} else {
  load(mwROOT("data","calculated","figure7b.RData"))
}

Plotting

Render and plot the figure.

plotdata <- bind_rows(plotdata, plotdatab)
levels(plotdata$type) <- c("Exponential", "Linear")
plotdata$type <- factor(plotdata$type, levels = c("Linear", "Exponential"))
sitedata$site <- factor(sitedata$site, levels = c(sitenames[-1], "Combined"))

# ggplot color palette
gg_color_hue <- function(n) {
    hues = seq(15, 375, length = n + 1)
    hcl(h = hues, l = 65, c = 100)[1:n]
}
gg_colors <- gg_color_hue(5)

plotdata %>% ggplot(aes(x = meanlog,
                        y = sdlog,
                        z = growth_rate)) +
  geom_raster(aes(fill = growth_rate)) +
  geom_contour(colour = "white", alpha = 0.8) + 
  scale_fill_gradientn("Growth\nRate", 
                           colours=c("#00000000","#BBBBBBBB"),
                           limits=c(min(plotdata$growth_rate), 
                                    max(plotdata$growth_rate))) +
  facet_grid(. ~ type) + 
  scale_x_continuous(limits = c(min(meanlog), 
                                max(meanlog)), 
                     expand = c(0, 0)) +
  scale_y_continuous(limits = c(min(sdlog), 
                                max(sdlog)), 
                     expand = c(0, 0)) +
  xlab(expression(paste("Log-normal location parameter (", mu, ")"))) + 
  ylab(expression(paste("Log-normal scale parameter (", sigma, ")"))) +
  labs(fill = "Growth\nRate", color = "Sites        ") +    # Weird bug requiring spaces to order legends 
  geom_point(aes(x = meanlog,
                 y = sdlog,
                 z = NA,
                 color = site), 
             size = 1.8,
             data = sitedata) +
  scale_color_manual(values = c(gg_colors, "#000000")) +
  theme_bw() ->
  p

myline <- lm(sdlog ~ meanlog, data = sitedata)
plotdata %>%
  ggplot(aes(x = meanlog,
             y = sdlog,
             z = growth_rate)) +
  geom_raster(aes(fill = growth_rate)) +
  geom_contour(colour = "white", alpha = 0.8) +
  geom_point(aes(color = site, z = NA), data=sitedata) +
  geom_abline(slope = myline$coefficients['meanlog'], 
              intercept = myline$coefficients['(Intercept)'],
              color = "red") ->
  p1
p1

Save figure.

ggsave("Figure7_LinearVsNonlinear.png", width=7, height=4, device = "png", p)

mdlama/milkweed documentation built on Sept. 15, 2022, 9:24 a.m.