In mdlama/milkweed:

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

Random effect error bar reference: https://stackoverflow.com/questions/13847936/in-r-plotting-random-effects-from-lmer-lme4-package-using-qqmath-or-dotplot

# Need to rescale standard errors!!!!
# For growth, sdr = attr(scaled$h_apical.next, 'scaled:scale'), else, = 1
# Let mua = attr(scaled$h_apical, 'scaled:center'), 
#     sda = attr(scaled$h_apical, 'scaled:scale'),
#     sdh = attr(scaled$herb_avg, 'scaled:scale')
# Then, seh_unscaled = seh_scaled*sdr/sdh - sea_scaled*sdr*mua/(sda*sdh).
# Note that sea_scaled is zero when there is no interaction term.
# Key: **h = herbivory; **a = apical height

years <- ipm$all_years
mydata <- tbl_df(data.frame())

Translated conditional modes

Survival

attach(ipm$pars)

# Get se from random effect
re <- lme4::ranef(surv.fit$mdl, condVar = TRUE, whichel = "year")[[1]]
pv <- attr(re, "postVar")
if (any(colnames(re) == "herb_avg")) {
  ind <- which(colnames(re) == "herb_avg")
  se <- sqrt(pv[ind,ind,])  # Index corresponds to slope parameter of interest
} else {
  se <- rep(0,length(years))
}

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

# Only necessary for adjustments based on interaction term
qsadj <- surv.fit$pars$unscaled["Bertha", "h_apical:herb_avg"]*quantile(h_apical, probs = c(0.05, 0.5, 0.95))

# Error bars for herb_avg random effect (NO interaction term)
sdr <- 1
mua <- attr(surv.fit$scaled$h_apical, "scaled:center")
sda <- attr(surv.fit$scaled$h_apical, "scaled:scale")
sdh <- attr(surv.fit$scaled$herb_avg, "scaled:scale")
ci <- 1.96*se*sdr/sdh

# Adjust error bars for interaction term if present
if (any(rownames(summary(surv.fit$mdl)$coefficients) == "h_apical:herb_avg")) {
  sei <- summary(surv.fit$mdl)$coefficients["h_apical:herb_avg", "Std. Error"]
  cis <- 1.96*sei*mua*sdr/(sda*sdh)
} else {
  cis <- 0;
}
ci <- ci - cis # Implement adjustment

mydata <- bind_rows(mydata, data.frame(rate = "Survival Probability",
                                       year = years, 
                                       slope = surv.fit$pars$unscaled[years,"herb_avg"] + qsadj["50%"],
                                       ci = ci))
detach(ipm$pars)

Flowering

attach(ipm$pars)

# Get se from random effect
re <- lme4::ranef(flower.fit$mdl, condVar = TRUE, whichel = "year")[[1]]
pv <- attr(re, "postVar")
if (any(colnames(re) == "herb_avg")) {
  ind <- which(colnames(re) == "herb_avg")
  se <- sqrt(pv[ind,ind,])  # Index corresponds to slope parameter of interest
} else {
  se <- rep(0,length(years))
}

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

# Only necessary for adjustments based on interaction term
qfadj <- flower.fit$pars$unscaled["Bertha", "h_apical:herb_avg"]*quantile(h_apical, probs = c(0.05, 0.5, 0.95))

# Error bars for herb_avg random effect (NO interaction term)
sdr <- 1
mua <- attr(flower.fit$scaled$h_apical, "scaled:center")
sda <- attr(flower.fit$scaled$h_apical, "scaled:scale")
sdh <- attr(flower.fit$scaled$herb_avg, "scaled:scale")
ci <- 1.96*se*sdr/sdh

# Adjust error bars for interaction term if present
if (any(rownames(summary(flower.fit$mdl)$coefficients) == "h_apical:herb_avg")) {
  sei <- summary(flower.fit$mdl)$coefficients["h_apical:herb_avg", "Std. Error"]
  cif <- 1.96*sei*mua*sdr/(sda*sdh)
} else {
  cif <- 0;
}
ci <- ci - cif # Implement adjustment

mydata <- bind_rows(mydata, data.frame(rate = "Flowering Probability",
                                       year = years,
                                       slope = flower.fit$pars$unscaled[years,"herb_avg"] + qfadj["50%"],
                                       ci = ci))
detach(ipm$pars)

Growth

attach(ipm$pars)

# Get se from random effect
re <- lme4::ranef(growth.fit$mdl, condVar = TRUE, whichel = "year")[[1]]
pv <- attr(re, "postVar")
if (any(colnames(re) == "herb_avg")) {
  ind <- which(colnames(re) == "herb_avg")
  se <- sqrt(pv[ind,ind,])  # Index corresponds to slope parameter of interest
} else {
  se <- rep(0,length(years))
}

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

# Only necessary for adjustments based on interaction term
qgadj <- growth.fit$pars$unscaled["Bertha", "h_apical:herb_avg"]*quantile(h_apical, probs = c(0.05, 0.5, 0.95))

# Error bars for herb_avg random effect (NO interaction term)
sdr <- attr(growth.fit$scaled$h_apical.next, "scaled:scale")
mua <- attr(growth.fit$scaled$h_apical, "scaled:center")
sda <- attr(growth.fit$scaled$h_apical, "scaled:scale")
sdh <- attr(growth.fit$scaled$herb_avg, "scaled:scale")
ci <- 1.96*se*sdr/sdh

# Adjust error bars for interaction term if present
if (any(rownames(summary(growth.fit$mdl)$coefficients) == "h_apical:herb_avg")) {
  sei <- summary(growth.fit$mdl)$coefficients["h_apical:herb_avg", "Std. Error"]
  cig <- 1.96*sei*mua*sdr/(sda*sdh)
} else {
  cig <- 0;
}
ci <- ci - cig # Implement adjustment

mydata <- bind_rows(mydata, data.frame(rate = "Growth",
                                       year = years, 
                                       slope = growth.fit$pars$unscaled[years,"herb_avg"] + qgadj["50%"],
                                       ci = ci))
detach(ipm$pars)

Pods

attach(ipm$pars)

# Get se from random effect
re <- lme4::ranef(pods.fit$mdl, condVar = TRUE, whichel = "year")[[1]]
pv <- attr(re, "postVar")
if (any(colnames(re) == "herb_avg")) {
  ind <- which(colnames(re) == "herb_avg")
  se <- sqrt(pv[ind,ind,])  # Index corresponds to slope parameter of interest
} else {
  se <- rep(0,length(years))
}

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

# Only necessary for adjustments based on interaction term
qpadj <- pods.fit$pars$unscaled["Bertha", "h_apical.next:herb_avg"]*quantile(h_apical.next, probs = c(0.05, 0.5, 0.95))

# Error bars for herb_avg random effect (NO interaction term)
sdr <- 1
mua <- attr(pods.fit$scaled$h_apical.next, "scaled:center")
sda <- attr(pods.fit$scaled$h_apical.next, "scaled:scale")
sdh <- attr(pods.fit$scaled$herb_avg, "scaled:scale")
ci <- 1.96*se*sdr/sdh

# Adjust error bars for interaction term if present
if (any(rownames(summary(pods.fit$mdl)$coefficients) == "h_apical.next:herb_avg")) {
  sei <- summary(pods.fit$mdl)$coefficients["h_apical.next:herb_avg", "Std. Error"]
  cip <- 1.96*sei/sdh
} else {
  cip <- 0
}

mydata <- bind_rows(mydata, data.frame(rate = "Pod Production",
                                       year = years, 
                                       slope = pods.fit$pars$unscaled[years,"herb_avg"] + qpadj["50%"],
                                       ci = ci))
detach(ipm$pars)

Fixed effect vertical lines

attach(ipm$pars)
rdata <- data.frame(rate = c("Survival Probability",
                             "Flowering Probability",
                             "Growth",
                             "Pod Production"),
                    rslope = c(surv.fit$pars$unscaled["Bertha", "herb_avg"]+qsadj["50%"],
                               flower.fit$pars$unscaled["Bertha", "herb_avg"]+qfadj["50%"],
                               growth.fit$pars$unscaled["Bertha", "herb_avg"]+qgadj["50%"],
                               pods.fit$pars$unscaled["Bertha", "herb_avg"])+qpadj["50%"])

# Fixed effect transparent bars
fdata <- data.frame(rate = c("Survival Probability",
                             "Flowering Probability",
                             "Growth",
                             "Pod Production"),
                    fslope = rdata$rslope,
                    ferr = c(cis, cif, cig, cip))

# Hidden points to set axes for facets (shouldn't need if using grid.arrange)
lim_data <- data.frame(rate = c(rep("Survival Probability", 2),
                                rep("Flowering Probability", 2),
                                rep("Growth", 2),
                                rep("Pod Production", 2)),
                       lims = c(-1, 1,
                                0+qfadj[2], 3.5+qfadj[2],
                                -3+qgadj[2], 15+qgadj[2],
                                -1, 0.2))
detach(ipm$pars)

Render and save plot

title.size <- 10

Flowering

type = "Flowering Probability"
pf <- mydata %>% filter(rate == type) %>%
           ggplot(aes(x = slope, 
                      y = year)) + 
           geom_point() +
           geom_errorbarh(aes(xmin = slope - ci, 
                              xmax = slope + ci), 
                          height=0) +
           geom_vline(data=filter(rdata, rate == type), 
                      aes(xintercept = rslope), 
                      linetype=2) + 
           scale_x_continuous(limits = c(-2, 2)) +
           theme_bw() +
           ggtitle(paste("(a)", type)) +
           theme(plot.title = element_text(family = "Trebuchet MS", 
                                           size = title.size,
                                           hjust = 0.5),
                 axis.title.y = element_blank(),
                 axis.title.x = element_blank())

if (any(qfadj > 0)) {
  pf <- pf + annotate("point", 
                      x = qfadj[1]-qfadj[2], 
                      y = 0.55, 
                      shape=25) +
    annotate("point", 
             x = qfadj[3]-qfadj[2], 
             y = 0.55, 
             shape=25, 
             fill="black")
}

pf

Growth

type = "Growth"
pg <- mydata %>% filter(rate == type) %>%
           ggplot(aes(x = slope, 
                      y = year)) + 
           geom_point() +
           geom_errorbarh(aes(xmin = slope - ci, 
                              xmax = slope + ci), 
                          height=0) +
           geom_vline(data=filter(rdata, rate == type), 
                      aes(xintercept = rslope), 
                      linetype = 2) + 
           scale_x_continuous(limits = c(-15, 15)) + 
           theme_bw() +
           ggtitle(paste("(b)", type)) +
           theme(plot.title = element_text(family = "Trebuchet MS", 
                                           size = title.size,
                                           hjust = 0.5),
                 axis.title.y = element_blank(),
                 axis.text.y = element_blank(),
                 axis.ticks.y = element_blank(),
                 axis.title.x = element_blank())

if (any(qgadj > 0)) {
  pg <- pg + annotate("point", 
                      x = qgadj[1]-qgadj[2], 
                      y = 0.55, 
                      shape=25) +
    annotate("point", 
             x = qgadj[3]-qgadj[2], 
             y = 0.55, 
             shape=25, 
             fill="black")
}

pg

Survival

type = "Survival Probability"
ps <- mydata %>% filter(rate == type) %>%
           ggplot(aes(x = slope, 
                      y = year)) + 
           geom_point() +
           geom_errorbarh(aes(xmin = slope - ci, 
                              xmax = slope + ci), 
                          height=0) +
           geom_vline(data=filter(rdata, rate == type), 
                      aes(xintercept = rslope), 
                      linetype = 2) + 
           scale_x_continuous(limits = c(-2, 2)) +
           theme_bw() +
           ggtitle(paste("(d)", type)) +
           theme(plot.title = element_text(family = "Trebuchet MS", 
                                           size = title.size,
                                           hjust = 0.5),
                 axis.title.y = element_blank(),
                 axis.text.y = element_blank(),
                 axis.ticks.y = element_blank(),
                 axis.title.x = element_blank())

if (any(qsadj > 0)) {
  ps <- ps + annotate("point", 
                      x = qsadj[1]-qsadj[2], 
                      y = 0.55, 
                      shape=25) +
    annotate("point", 
             x = qsadj[3]-qsadj[2], 
             y = 0.55, 
             shape=25, 
             fill="black")
}

ps

Pods

type = "Pod Production"
pp <- mydata %>% filter(rate == type) %>%
           ggplot(aes(x = slope, 
                      y = year)) + 
           geom_point() +
           geom_errorbarh(aes(xmin = slope - ci, 
                              xmax = slope + ci), 
                          height = 0) +
           geom_vline(data=filter(rdata, rate == type), 
                      aes(xintercept = rslope), 
                      linetype = 2) + 
           scale_x_continuous(limits = c(-2.0, 2.0)) +
           theme_bw() +
           ggtitle(paste("(c)", type)) +
           theme(plot.title = element_text(family = "Trebuchet MS", 
                                           size = title.size,
                                           hjust = 0.5),
                 axis.title.y = element_blank(),
                 axis.title.x = element_blank())

if (any(qpadj > 0)) {
  pp <- pp + annotate("point", 
                      x = qpadj[1]-qpadj[2], 
                      y = 0.55, 
                      shape=25) +
    annotate("point", 
             x = qpadj[3]-qpadj[2], 
             y = 0.55, 
             shape=25, 
             fill="black")
}

pp

All together now

pall <- gridExtra::grid.arrange(pf, pg, pp, ps, 
                     left="Year", 
                     bottom=grid::textGrob(latex2exp::TeX('$\\beta_{z_{\\omega}}$')))
ggsave("Figure3_HerbSlopePars.png", width=6, device = "png", pall)