In mdlama/milkweed:

Initialization

library(milkweed)
library(dplyr)
library(ggplot2)
library(lme4)
library(AICcmodavg)
requirePackages("gridExtra")

stemdata$year <- factor(stemdata$year)

Fit a random effects model on probability of getting eaten ($\mathrm{p}_{\omega}$)

First, let's get the data ready. We need to normalize the data (mean 0 and standard deviation 1).

stemdata_usc <- stemdata %>% filter(!is.na(h_apical),
                                    !is.na(herb_avg))


stemdata_sc <- stemdata_usc %>% mutate_at(.vars = vars(h_apical, herb_avg),
                                          .funs = funs(as.numeric(scale(.))))

# Get scaling parameters (really just mean and sd)
h_apical_sc <- scale(stemdata_usc$h_apical)
herb_avg_sc <- scale(stemdata_usc$herb_avg)

glmerCtrl <- glmerControl(optimizer = c("bobyqa"), optCtrl = list(maxfun=50000))

Let's choose random effects first - we'll focus on intercept and slope on apical height, ignoring herbivory slopes.

# Random effects first
herbivory.reg.mdl1 <- glmer(munched ~ h_apical + (h_apical|site/transect) + (h_apical|year), data=stemdata_sc, nAGQ=1, family=binomial(), control = glmerCtrl)

herbivory.reg.mdl2 <- glmer(munched ~ h_apical + (1|site/transect) + (h_apical|year), data=stemdata_sc, nAGQ=1, family=binomial(), control = glmerCtrl)

herbivory.reg.mdl3 <- glmer(munched ~ h_apical + (h_apical|year), data=stemdata_sc, nAGQ=1, family=binomial(), control = glmerCtrl)

herbivory.reg.mdl4 <- glmer(munched ~ h_apical + (1|year), data=stemdata_sc, nAGQ=1, family=binomial(), control = glmerCtrl)

herbivory.reg.mdl5 <- glmer(munched ~ h_apical + (1|site/transect) + (1|year), data=stemdata_sc, nAGQ=1, family=binomial(), control = glmerCtrl)

aictab(list(herbivory.reg.mdl1, 
            herbivory.reg.mdl2, 
            herbivory.reg.mdl3,
            herbivory.reg.mdl4,
            herbivory.reg.mdl5))

bictab(list(herbivory.reg.mdl1, 
            herbivory.reg.mdl2, 
            herbivory.reg.mdl3,
            herbivory.reg.mdl4,
            herbivory.reg.mdl5))

herbivory.reg.winner <- herbivory.reg.mdl1

Winner is model with site as a random effect on intercept only, and year as a random effect on intercept and apical height (going with BIC).

munched ~ h_apical + (h_apical|site/transect) + (h_apical|year)

Let's check the fixed effects. In this case, we just want to compare Model #1 to the model with no apical height.

herbivory.reg.noapical <- glmer(munched ~ 1 + (h_apical|site/transect) + (h_apical|year), data=stemdata_sc, nAGQ=1, family=binomial(), control = glmerCtrl)

aictab(list(herbivory.reg.winner, 
            herbivory.reg.noapical),
       modnames=c("Height", "No Height"))

bictab(list(herbivory.reg.winner, 
            herbivory.reg.noapical),
       modnames=c("Height", "No Height"))

As the model with no apical height is comparable, we use the simpler model. Thus, we move forward under the assumption that apical height does not predict probability of being eaten.

Let's look at a plot to help us see the variation across site and year.

binomial_smooth <- function(...) {
  geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}

pa <- stemdata_usc %>% ggplot(aes(x = h_apical, 
                                  y = munched)) + 
  geom_jitter(alpha = 0.3,
              height = 0.1) +
  facet_grid(year ~ site) +
  binomial_smooth() +
  theme_bw() +
  xlab("Height (cm)") +
  ylab("Probability of being eaten") +
  ggtitle("(a)") + 
  theme(plot.title = element_text(hjust = 0,
                                  margin = margin(b = 15, unit = "pt")))

pa

Fit a random effects model on herbivory magnitude

Same approach as above, i.e. let's select random effects first. This time, going by Zuur, we need to specify REML=T. Also note that we need to condition the data to only those stems that experience herbivory (filter(munched == 1)).

herbivory.reg.mdl1 <- lmer(herb_avg ~ h_apical + (h_apical|site/transect) + (h_apical|year), data=stemdata_sc %>% filter(munched == 1), REML=T)

herbivory.reg.mdl2 <- lmer(herb_avg ~ h_apical + (1|site/transect) + (h_apical|year), data=stemdata_sc %>% filter(munched == 1), REML=T)

herbivory.reg.mdl3 <- lmer(herb_avg ~ h_apical + (h_apical|site/transect) + (1|year), data=stemdata_sc %>% filter(munched == 1), REML=T)

herbivory.reg.mdl4 <- lmer(herb_avg ~ h_apical + (h_apical|year), data=stemdata_sc %>% filter(munched == 1), REML=T)

herbivory.reg.mdl5 <- lmer(herb_avg ~ h_apical + (h_apical|site/transect), data=stemdata_sc %>% filter(munched == 1), REML=T)

aictab(list(herbivory.reg.mdl1, 
            herbivory.reg.mdl2,
            herbivory.reg.mdl3,
            herbivory.reg.mdl4,
            herbivory.reg.mdl5))

bictab(list(herbivory.reg.mdl1, 
            herbivory.reg.mdl2,
            herbivory.reg.mdl3,
            herbivory.reg.mdl4,
            herbivory.reg.mdl5))

Looks like the "full" model for random effects. Now on to fixed effects. For the model selection part, we set REML=F. Note the need to recompute for the random effects winner.

herbivory.reg.winner <- lmer(herb_avg ~ h_apical + (h_apical|site/transect) + (h_apical|year), data=stemdata_sc %>% filter(munched == 1), REML=F)

herbivory.reg.noapical <- lmer(herb_avg ~ 1 + (h_apical|site/transect) + (h_apical|year), data=stemdata_sc %>% filter(munched == 1), REML=F)

aictab(list(herbivory.reg.winner, 
            herbivory.reg.noapical),
       modnames = c("Height", "No Height"))

bictab(list(herbivory.reg.winner, 
            herbivory.reg.noapical),
       modnames = c("Height", "No Height"))

Similar to probability of being eaten, it appears that the model with no apical height is comparable to the model with apical height. For parsimony, we will assume that apical height does not predict magnitude of herbivory.

Let's look at a plot for some visual intuition.

pb <- stemdata_usc %>% 
  filter(munched == 1) %>%
  ggplot(aes(x = h_apical, 
             y = herb_avg)) + 
  facet_grid(year ~ site) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm") +
  theme_bw() +
  xlab("Height (cm)") +
  ylab("Herbivory score") +
  ggtitle("(b)") + 
  theme(plot.title = element_text(hjust = 0,
                                  margin = margin(b = 15, unit = "pt")))

pb

Now we combine both plots and save.

p <- gridExtra::grid.arrange(pa, pb, ncol=1)
ggsave("AppendixFigure3_HerbVsHeight.png", width=6, height=8, device = "png", p)

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