In mdlama/milkweed:

Load necessary packages and data and perform initializations.

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

data("stemdata")
metadata <- stemdata
metadata$year <- factor(metadata$year)
metadata$transect <- factor(metadata$transect)

# https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q4/022870.html
glmerCtrl <- glmerControl(optimizer = c("bobyqa"), optCtrl = list(maxfun=50000))
lmerCtrl <- lmerControl(optimizer = c("bobyqa"), optCtrl = list(maxfun=50000))

Model Selection

We will use the top-down approach for model selection discussed in Zuur et al 2009 (reccomended by Diggle et al 2002) starting with the full ("beyond optimal") model and refining the random effects first before moving on to the fixed effects.

Flowering

metadata_usc <- metadata %>% filter(!is.na(h_apical),
                                    !is.na(herb_avg),
                                    !is.na(fec.flower))

metadata_sc <- metadata_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(metadata_usc$h_apical)
herb_avg_sc <- scale(metadata_usc$herb_avg)

metadata_sc %>% ggplot(aes(x = h_apical, 
                           y = fec.flower, 
                           color = year)) + 
                geom_jitter(height = 0.1, 
                            alpha = 0.3) +
                facet_wrap(~ site)

Random Effects

Use nAGQ = 0 for quick approximations, but less exact form of parameter estimation.

flower.reg.mdl1 <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical*herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

Let's try and simplify site/transect random effect first, but we'll keep herb_avg in as that is our main interest.

flower.reg.mdl2=glmer(fec.flower ~ h_apical*herb_avg + (herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

flower.reg.mdl3=glmer(fec.flower ~ h_apical*herb_avg + (1|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

aictab(list(flower.reg.mdl1, flower.reg.mdl2, flower.reg.mdl3), 
       modnames = c("mdl1", "mdl2", "mdl3"))

Models 2 and 3 are indistinguishable based on AIC, so we will move onto fixed effect selection keeping the more complicated random effect model, leading us to choose model 2. Let's move to year and compare all the different models.

flower.reg.mdl4=glmer(fec.flower ~ h_apical*herb_avg + (herb_avg|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())
flower.reg.mdl5=glmer(fec.flower ~ h_apical*herb_avg + (herb_avg|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=0, family=binomial())
flower.reg.mdl6=glmer(fec.flower ~ h_apical*herb_avg + (herb_avg|site/transect) + (herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())
flower.reg.mdl7=glmer(fec.flower ~ h_apical*herb_avg + (herb_avg|site/transect) + (1|year), data=metadata_sc, nAGQ=0, family=binomial())

aictab(list(flower.reg.mdl2, flower.reg.mdl4, flower.reg.mdl5, flower.reg.mdl6, flower.reg.mdl7),
       modnames = c("mdl2","mdl4","mdl5","mdl6","mdl7"))

# AIC still indicates Models 2 easily

# nAGQ=1 uses Laplace approximation rather than quick approximation of before, for better estimation in the winning model.
flower.reg.winner=glmer(fec.flower ~ h_apical*herb_avg + (herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

summary(flower.reg.winner)

Fixed Effects

Moving onto selection for fixed effects, first drop h_apical:herb_avg interaction term, since the p-value is the lowest non-stat. sig. then eliminate lowest term by term until only stat. sig. terms remain.

flower.reg.mdl2a = glmer(fec.flower ~ h_apical + herb_avg + (herb_avg|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)
summary(flower.reg.mdl2a)

flower.reg.mdl2b = glmer(fec.flower ~ h_apical + herb_avg - 1 + (herb_avg|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)
summary(flower.reg.mdl2b)

flower.reg.mdl2c = glmer(fec.flower ~ h_apical - 1 + (herb_avg|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)
summary(flower.reg.mdl2c)

aictab(list(flower.reg.winner, flower.reg.mdl2a, flower.reg.mdl2b, flower.reg.mdl2c), 
       modnames = c("winner", "mdl2a", "mdl2b", "mdl2c"))

Looks like Model 2c and 2b are indistinguishable based on AICc. In the event of ties in the fixed effects, we have decided to use the model that includes herbivory, since we are interested to see if it is a predictor. If it is only a weak predictor after all, we will see that in sensitivity analysis of the IPM. Let's check p-values and create class to verify that unscaled version is correct.

flower.reg.final <- flower.reg.mdl2b
summary(flower.reg.final)

Model Validation and Diagnostics

Now some diagnostics (aka model validation). Creating the standard residual plots. Interpretting these plots is questionable for a non-linear model fit, so we will also imploy non-linear diagnostics next.

par(mfrow=c(2,2))
plot(fitted(flower.reg.final), residuals(flower.reg.final, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(flower.reg.final, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(flower.reg.final, "deviance"), xlab="Herbivory average", ylab="Deviance residuals")
par(mfrow=c(1,1))

DHARMa package diagnostics for non-linear, hierarchical, mixed models. This simulates new data from the fitted model for the predictor variable combination of each observation. Then for each observation, calculates the empirical cumulative density function for the simulated data, which describes the expected spread for an observation at the respective point in predictor space, conditional on the fitted model. The residual is defined as the value of the empirical density function at the value of the observed data (DHARMa package vignette).

flower.simulationOutput <- DHARMa::simulateResiduals(fittedModel = flower.reg.final, n=1000) #simulate data, n set to 1000 for high precision

# DHARMa diagnostic plot, then plots of residuals against the predictor variables independently, we want flat horizontal lines at 0.25, 0.5, and 0.75 in a perfect world.
DHARMa::plotSimulatedResiduals(flower.simulationOutput)
DHARMa::plotResiduals(metadata_sc$h_apical, flower.simulationOutput$scaledResiduals)
DHARMa::plotResiduals(metadata_sc$herb_avg, flower.simulationOutput$scaledResiduals)

# tests to see if the residuals fit the expected distribution
DHARMa::testUniformity(flower.simulationOutput)

Plots look great! Looks like the final model is a good fit. The Uniformity Test is probably failing since we have such a large sample size and it is not a particularly strong test. It does have a relatively large p-value compared to some other decent fits we've seen.

Visualizations

Create model class

scaled <- list(h_apical = h_apical_sc, herb_avg = herb_avg_sc)
flower.fit <- mwMod(list(mdl = flower.reg.final, vars = c("h_apical", "herb_avg"), scaled = scaled))
checkPars(flower.fit) # Check parameters

Scaled parameters across site:

flower.fit$pars$scaled

Unscaled parameters across site:

flower.fit$pars$unscaled

Let's plot full data.

# Add points
flower.plot <- metadata_usc %>% 
  ggplot(aes(x = h_apical, 
             y = fec.flower, 
             color = site, 
             size = herb_avg)) +
  geom_jitter(width = 0, 
              height = 0.05,
              alpha = 0.4)

# Create plot data for prediction curves
plotdata <- metadata_usc %>% 
  group_by(site) %>% 
  summarize(mean_herb_avg = mean(herb_avg),
            min_h_apical = min(h_apical),
            max_h_apical = max(h_apical)) %>%
  add_row(site = "Bertha",
          mean_herb_avg = mean(metadata_usc$herb_avg, na.rm=T),
          min_h_apical = min(metadata_usc$h_apical, na.rm=T),
          max_h_apical = max(metadata_usc$h_apical, na.rm=T))

# Create prediction curves
N <- 100
suppressWarnings( # Hiding the follow warnings: (1) Unequal factor levels: coercing to character and (2) binding character and factor vector, coercing into character vector
  curves <- bind_rows(
    lapply(as.vector(plotdata$site), function(site) {
      i <- which(plotdata$site == site)
      data.frame(
        site = site,
        h_apical = seq(plotdata$min_h_apical[i], plotdata$max_h_apical[i], length.out = N),
        herb_avg = plotdata$mean_herb_avg[i],
        fec.flower = predict(
          flower.fit,
          type = site,
          newdata = data.frame(
            h_apical = seq(plotdata$min_h_apical[i], plotdata$max_h_apical[i], length.out = N),
            herb_avg = rep(plotdata$mean_herb_avg[i], N)
          )
        )
      )
    })
  )
)

flower.plot + geom_line(data = curves)

Now site-specific models.

flower.plot <- metadata_usc %>% 
  ggplot(aes(x = h_apical, 
             y = fec.flower)) + 
  geom_jitter(height = 0.1, 
              alpha = 0.3) +
  facet_wrap(~ site)

flower.plot + 
  geom_line(data = curves %>% filter(site != "Bertha")) + 
  facet_wrap(~ site)

Survival

metadata_usc <- metadata %>% filter(!is.na(h_apical),
                                    !is.na(herb_avg),
                                    fec.flower == 1, # filtering out 0s because only the sexual pathway uses the survival function
                                    !is.na(surv))

metadata_sc <- metadata_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(metadata_usc$h_apical)
herb_avg_sc <- scale(metadata_usc$herb_avg)

metadata_sc %>% ggplot(aes(x = h_apical, 
                           y = surv, 
                           color = year)) + 
                geom_jitter(height = 0.1, 
                            alpha = 0.3) +
                facet_wrap(~ site)

Random Effects

First fit full model. Using nAGQ = 0 for quick approximation (matches results with nAGQ = 1).

surv.reg.full = glmer(surv ~ h_apical*herb_avg +  (h_apical*herb_avg|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, nAGQ = 0, family=binomial())
summary(surv.reg.full)

First thing to go is interaction on site and transect, since it has very low variance and is contributing relatively little to the variation in the model.

surv.reg.mdl1 = glmer(surv ~ h_apical*herb_avg +  (h_apical+herb_avg|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())
summary(surv.reg.mdl1)

Interaction on year is also minimal in both, so it will be removed next.

surv.reg.mdl2 = glmer(surv ~ h_apical*herb_avg +  (h_apical+herb_avg|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())
summary(surv.reg.mdl2)

Dropping herb_avg in site/transect next.

surv.reg.mdl3 = glmer(surv ~ h_apical*herb_avg +  (h_apical | site/transect) + (h_apical + herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())
summary(surv.reg.mdl3)

Now drop herb_avg on year.

surv.reg.mdl4 = glmer(surv ~ h_apical*herb_avg +  (h_apical | site/transect) + (h_apical | year), data=metadata_sc, nAGQ=0, family=binomial())
summary(surv.reg.mdl4)

Now all that's left is to try combos of dropping h_apical.

surv.reg.mdl5 = glmer(surv ~ h_apical*herb_avg +  (1|site/transect)+(h_apical|year), data=metadata_sc, nAGQ=0, family=binomial())
surv.reg.mdl6 = glmer(surv ~ h_apical*herb_avg +  (h_apical|site/transect)+(1|year), data=metadata_sc, nAGQ=0, family=binomial())
surv.reg.mdl7 = glmer(surv ~ h_apical*herb_avg +  (1|site/transect)+(1|year), data=metadata_sc, nAGQ=0, family=binomial())

AIC table.

aictab(list(surv.reg.full, surv.reg.mdl1, surv.reg.mdl2, surv.reg.mdl3, surv.reg.mdl4, surv.reg.mdl5, surv.reg.mdl6, surv.reg.mdl7), modnames = c('surv.reg.full', 'surv.reg.mdl1', 'surv.reg.mdl2', 'surv.reg.mdl3', 'surv.reg.mdl4', 'surv.reg.mdl5', 'surv.reg.mdl6', 'surv.reg.mdl7'))

Winner: nearly indistinguishable AICc between Models 2 and 3. Let's go with Model 3 as there appears to be convergence issues with Model 2.

surv.reg.winner = glmer(surv ~ h_apical*herb_avg +  (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)
summary(surv.reg.winner)

Fixed Effects

On to fixed effects, removing h_apical first.

surv.reg.mdl1 = glmer(surv ~ herb_avg + h_apical:herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, family=binomial(), nAGQ=1, control=glmerCtrl)

surv.reg.mdl2 = glmer(surv ~ herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, family=binomial(), nAGQ=1, control=glmerCtrl)

surv.reg.mdl3 = glmer(surv ~ 1 + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, family=binomial(), nAGQ=1, control=glmerCtrl)

sAICtab <- AIC(surv.reg.winner, surv.reg.mdl1, surv.reg.mdl2, surv.reg.mdl3)
sAICtab$dAIC <- min(sAICtab$AIC)-sAICtab$AIC
sAICtab

sBICtab <- BIC(surv.reg.winner, surv.reg.mdl1, surv.reg.mdl2, surv.reg.mdl3)
sBICtab$dBIC <- min(sBICtab$BIC)-sBICtab$BIC
sBICtab

# mdl3 wins by large margin in BIC but close in AIC to mdl2 - since we're interested in herbivory, let's go with mdl2
surv.reg.final <- surv.reg.mdl2
summary(surv.reg.final)

Model Validation and Diagnostics

Now some diagnostics (aka model validation). Creating standard residual plots.

par(mfrow=c(2,2))
plot(fitted(surv.reg.final), residuals(surv.reg.final, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(surv.reg.final, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(surv.reg.final, "deviance"), xlab="Herbivory average", ylab="Deviance residuals")
par(mfrow=c(1,1))

DHARMa Diagnostics. Same process as above.

surv.simulationOutput = DHARMa::simulateResiduals(fittedModel = surv.reg.final, n=1000)

DHARMa::plotSimulatedResiduals(surv.simulationOutput)
DHARMa::plotResiduals(metadata_sc$h_apical, surv.simulationOutput$scaledResiduals)
DHARMa::plotResiduals(metadata_sc$herb_avg, surv.simulationOutput$scaledResiduals)

Here we actually pass the Uniformity test, indicating that the residuals do fit the expected distribution.

DHARMa::testUniformity(surv.simulationOutput)

# Create model class
scaled <- list(h_apical = h_apical_sc, herb_avg = herb_avg_sc)
surv.fit <- mwMod(list(mdl = surv.reg.final, vars = c("h_apical", "herb_avg"), scaled = scaled))
checkPars(surv.fit) # Check parameters

Growth

# Transform data
metadata_usc <- metadata %>% filter(!is.na(h_apical),
                                    !is.na(h_apical.next),
                                    !is.na(herb_avg),
                                    fec.flower == 1, #growth is also only used in the sexual pathway
                                    surv == 1) %>%   #plants must survive until september to reprodcude in our model
  mutate(log_herb_avg = log(0.01 + herb_avg), arcsin_herb_avg = asin(sqrt(herb_avg/6)), h_apical_sq = h_apical*h_apical)

metadata_sc <- metadata_usc %>% mutate_at(.vars = vars(h_apical, h_apical.next, herb_avg, log_herb_avg, arcsin_herb_avg, h_apical_sq),
                                          .funs = funs(as.numeric(scale(.))))

# Get scaling parameters (really just mean and sd)
h_apical_sc <- scale(metadata_usc$h_apical)
h_apical.next_sc <- scale(metadata_usc$h_apical.next)
herb_avg_sc <- scale(metadata_usc$herb_avg)
log_herb_avg_sc = scale(metadata_usc$log_herb_avg)
arcsin_herb_avg_sc = scale(metadata_usc$arcsin_herb_avg)
h_apical_sq_sc = scale(metadata_usc$h_apical_sq)

metadata_sc %>% ggplot(aes(x = h_apical, 
                           y = h_apical.next, 
                           color = year)) + 
                geom_point(alpha = 0.3) +
                facet_wrap(~ site)

metadata_sc %>% ggplot(aes(x = h_apical,
                           y = h_apical.next,
                           color = site)) + 
  geom_point(alpha = 0.3)

As this is a bona-fide LME, we can do the method mentioned in Zuur. First, model selection on random effects with REML.

growth.reg.full=lmer(h_apical.next~h_apical*herb_avg + (h_apical*herb_avg|site/transect)+(h_apical*herb_avg|year),data=metadata_sc, REML=T, control=lmerCtrl)
summary(growth.reg.full)

growth.reg.mdl1=lmer(h_apical.next~h_apical*herb_avg + (h_apical+herb_avg|site/transect)+(h_apical*herb_avg|year),data=metadata_sc, REML=T, control=lmerCtrl)
summary(growth.reg.mdl1)

growth.reg.mdl2=lmer(h_apical.next~h_apical*herb_avg + (h_apical|site/transect)+(h_apical*herb_avg|year),data=metadata_sc, REML=T, control=lmerCtrl)
summary(growth.reg.mdl2)

growth.reg.mdl3=lmer(h_apical.next~h_apical*herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year),data=metadata_sc, REML=T, control=lmerCtrl)
summary(growth.reg.mdl3)

growth.reg.mdl4=lmer(h_apical.next~h_apical*herb_avg + (1|site/transect)+(h_apical+herb_avg|year),data=metadata_sc, REML=T, control=lmerCtrl)
# growth.reg.mdl4=lmer(h_apical.next~h_apical*herb_avg + (h_apical|site/transect)+(herb_avg|year),data=metadata_sc, REML=T)
summary(growth.reg.mdl4)

growth.reg.mdl5=lmer(h_apical.next~h_apical*herb_avg + (1|site/transect)+(h_apical|year), data=metadata_sc, REML=T, control=lmerCtrl)
#growth.reg.mdl5=lmer(h_apical.next~h_apical*herb_avg + (h_apical|site/transect)+(1|year),data=metadata_sc, REML=T)
summary(growth.reg.mdl5)

growth.reg.mdl6=lmer(h_apical.next~h_apical*herb_avg + (1|site/transect)+(herb_avg|year), data=metadata_sc, REML=T, control=lmerCtrl)
#growth.reg.mdl6=lmer(h_apical.next~h_apical*herb_avg + (1|site/transect)+(1|year),data=metadata_sc, REML=T)
summary(growth.reg.mdl6)

aictab(list(growth.reg.full, growth.reg.mdl1, growth.reg.mdl2, growth.reg.mdl3, growth.reg.mdl4, growth.reg.mdl5, growth.reg.mdl6), modnames = c('growth.reg.full', 'growth.reg.mdl1', 'growth.reg.mdl2', 'growth.reg.mdl3', 'growth.reg.mdl4', 'growth.reg.mdl5', 'growth.reg.mdl6'))

# Model 2 wins AIC

Now fixed effects, with ML this time. Comparing all alternative models.

growth.reg.mdl2a=lmer(h_apical.next ~ h_apical*herb_avg + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)
#do single-term deletions to see which term needs to be dropped first
drop1(growth.reg.mdl2a, scope=c("h_apical", "herb_avg", "h_apical:herb_avg"))

# Looks like interaction goes first based on the drop1() results
growth.reg.mdl2b=lmer(h_apical.next ~ h_apical+herb_avg + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl) #all 3 remaining terms

growth.reg.mdl2c=lmer(h_apical.next ~ h_apical + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl) #2 remaining terms
growth.reg.mdl2d=lmer(h_apical.next ~ h_apical + herb_avg - 1 + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)
growth.reg.mdl2e=lmer(h_apical.next ~ herb_avg + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

growth.reg.mdl2f=lmer(h_apical.next ~ h_apical - 1 + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl) #1 remaining term
growth.reg.mdl2g=lmer(h_apical.next ~ herb_avg - 1 + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)
growth.reg.mdl2h=lmer(h_apical.next ~ 1 + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

aictab(list(growth.reg.mdl2a, growth.reg.mdl2b, growth.reg.mdl2c, growth.reg.mdl2d, growth.reg.mdl2e, growth.reg.mdl2f, growth.reg.mdl2g, growth.reg.mdl2h), modnames = c('growth.reg.mdl2a', 'growth.reg.mdl2b', 'growth.reg.mdl2c', 'growth.reg.mdl2d', 'growth.reg.mdl2e', 'growth.reg.mdl2f', 'growth.reg.mdl2g', 'growth.reg.mdl2h'))

# Models c, and b, are indistinguishable for AICc, so we will go with 2b.

Let's try with Model 1, as Model 2 above gives negative Sept height.

growth.reg.mdl1a=lmer(h_apical.next ~ h_apical*herb_avg + (h_apical+herb_avg|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)
#do single-term deletions to see which term needs to be dropped first
drop1(growth.reg.mdl1a, scope=c("h_apical", "herb_avg", "h_apical:herb_avg"))

# Looks like interaction goes first based on the drop1() results
growth.reg.mdl1b=lmer(h_apical.next ~ h_apical+herb_avg + (h_apical+herb_avg|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl) #all 3 remaining terms

growth.reg.mdl1c=lmer(h_apical.next ~ h_apical + (h_apical|site/transect)+(h_apical + herb_avg|year), data=metadata_sc, REML=F) #2 remaining terms
growth.reg.mdl1d=lmer(h_apical.next ~ h_apical + herb_avg - 1 + (h_apical|site/transect)+(h_apical + herb_avg|year), data=metadata_sc, REML=F)
growth.reg.mdl1e=lmer(h_apical.next ~ herb_avg + (h_apical|site/transect)+(h_apical + herb_avg|year), data=metadata_sc, REML=F)

growth.reg.mdl2f=lmer(h_apical.next ~ h_apical - 1 + (h_apical|site/transect)+(h_apical + herb_avg|year), data=metadata_sc, REML=F) #1 remaining term
growth.reg.mdl2g=lmer(h_apical.next ~ herb_avg - 1 + (h_apical|site/transect)+(h_apical + herb_avg|year), data=metadata_sc, REML=F)
growth.reg.mdl2h=lmer(h_apical.next ~ 1 + (h_apical|site/transect)+(h_apical + herb_avg|year), data=metadata_sc, REML=F)


aictab(list(growth.reg.mdl2a, growth.reg.mdl2b, growth.reg.mdl2c, growth.reg.mdl2d, growth.reg.mdl2e, growth.reg.mdl2f, growth.reg.mdl2g, growth.reg.mdl2h), modnames = c('growth.reg.mdl2a', 'growth.reg.mdl2b', 'growth.reg.mdl2c', 'growth.reg.mdl2d', 'growth.reg.mdl2e', 'growth.reg.mdl2f', 'growth.reg.mdl2g', 'growth.reg.mdl2h'))

# Models a, and b, are indistinguishable for AICc, so we will go with 2a.

Need to refit with REML.

# growth.reg.final = lmer(h_apical.next ~ h_apical*herb_avg + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=T)
growth.reg.final = lmer(h_apical.next ~ h_apical + (h_apical|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, REML=T)
summary(growth.reg.final)

Model Validation and Diagnostics

Now some diagnostics (aka model validation).

par(mfrow=c(2,2))
plot(fitted(growth.reg.final), residuals(growth.reg.final, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(growth.reg.final, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(growth.reg.final, "deviance"), xlab="Herbivory average", ylab="Deviance residuals")
par(mfrow=c(1,1))

DHARMa Diagnostics

growth.simulationOutput = DHARMa::simulateResiduals(fittedModel = growth.reg.final, n=1000)

DHARMa::plotSimulatedResiduals(growth.simulationOutput)
DHARMa::plotResiduals(metadata_sc$h_apical, growth.simulationOutput$scaledResiduals)
DHARMa::plotResiduals(metadata_sc$herb_avg, growth.simulationOutput$scaledResiduals)

Low p-value again may be due to large sample. Plots look ok if you exclude the tails and concentrate on where most of the data is centered, however the shape could be improved.

DHARMa::testUniformity(growth.simulationOutput)

Herb_avg transform

Let's try transforming herb_avg and putting that into the existing model structure. Going with log and an arcsine transform on the herbivory average scaled to be 0<herb_avg<1.

growth.reg.final.log = lmer(h_apical.next ~ h_apical*log_herb_avg + (h_apical|site/transect)+(h_apical*log_herb_avg|year), data=metadata_sc, REML=T)
summary(growth.reg.final.log)
growth.reg.final.arcsin = lmer(h_apical.next ~ h_apical*arcsin_herb_avg + (h_apical|site/transect)+(h_apical*arcsin_herb_avg|year), data=metadata_sc, REML=T)
summary(growth.reg.final.arcsin)

#both transformations produce a much lower AIC, with log coming out on top
aictab(list(growth.reg.final, growth.reg.final.log, growth.reg.final.arcsin), modnames = c('growth.reg.final', 'growth.reg.final.log', 'growth.reg.final.arcsin'))

Both have MUCH lower AICc than the final model, so let's see what model validation says.

par(mfrow=c(2,2))
plot(fitted(growth.reg.final.log), residuals(growth.reg.final.log, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(growth.reg.final.log, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(growth.reg.final.log, "deviance"), xlab="log(Herbivory average)", ylab="Deviance residuals")
par(mfrow=c(1,1))

par(mfrow=c(2,2))
plot(fitted(growth.reg.final.arcsin), residuals(growth.reg.final.arcsin, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(growth.reg.final.arcsin, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(growth.reg.final.arcsin, "deviance"), xlab="arcsin(sqrt(Herbivory average))", ylab="Deviance residuals")
par(mfrow=c(1,1))

None of these residual plots look any better, so we will see what DHARMa says.

growth.simulationOutput.log = DHARMa::simulateResiduals(fittedModel = growth.reg.final.log, n=1000)
growth.simulationOutput.arcsin = DHARMa::simulateResiduals(fittedModel = growth.reg.final.arcsin, n=1000)

DHARMa::plotSimulatedResiduals(growth.simulationOutput.log)
DHARMa::plotResiduals(metadata_sc$h_apical, growth.simulationOutput.log$scaledResiduals)
DHARMa::plotResiduals(metadata_sc$log_herb_avg, growth.simulationOutput.log$scaledResiduals)

DHARMa::plotSimulatedResiduals(growth.simulationOutput.arcsin)
DHARMa::plotResiduals(metadata_sc$h_apical, growth.simulationOutput.arcsin$scaledResiduals)
DHARMa::plotResiduals(metadata_sc$arcsin_herb_avg, growth.simulationOutput.arcsin$scaledResiduals)

These plots look marginally better at best. Also the QQ plots are no different. Given that the DHARMa diagnostics are at least similar, the only benfit we would get from transforming herb_avg would be a lower AICc, which would be questionable, given that there seem to be few other improvements. The aictab() function documentation warns against comparing model AICs between those that have transformed and untransformed response variables. While we do not transform the response variable, we believe that for the same reasons, it would not be ideal to compare models that have different transformations on the explanatories. Therefore, we will use our original final model as our best approximation, despite imperfect model validation results and a higher AICc. It should also be noted that we tried including a quadratic term for herb_avg, with similar results. The model selection for that as well a site by site and year by year break down validating the use of a linear model for growth can be found in the "Growth_Exploration.Rmd" script.

# Create model class
scaled <- list(h_apical = h_apical_sc, herb_avg = herb_avg_sc, h_apical.next = h_apical.next_sc)
growth.fit <- mwMod(list(mdl = growth.reg.final, vars = c("h_apical", "herb_avg", "h_apical.next"), scaled = scaled))
checkPars(growth.fit) # Check parameters

Pods

Data

metadata_usc <- metadata %>% filter(!is.na(h_apical),
                                    !is.na(h_apical.next),
                                    !is.na(herb_avg),
                                    fec.flower == 1, #pods only result from the sexual pathway
                                    surv == 1, #again, plants must survive to September
                                    !is.na(N_pods),
                                    h_apical.next > 50,
                                    !(site == "YTB" & h_apical <= 50)) #necessary to compensate for an issue with some of the plants in YTB site

metadata_sc <- metadata_usc %>% mutate_at(.vars = vars(h_apical, h_apical.next, herb_avg),
                                          .funs = funs(as.numeric(scale(.))))

# Get scaling parameters (really just mean and sd)
h_apical_sc <- scale(metadata_usc$h_apical)
h_apical.next_sc <- scale(metadata_usc$h_apical.next)
herb_avg_sc <- scale(metadata_usc$herb_avg)

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

metadata_sc %>% ggplot(aes(x = h_apical.next, 
                           y = N_pods, 
                           color = year)) + 
                geom_point(alpha = 0.3) +
                facet_wrap(~ site)

June vs. September Height as a Predictor

First, do full model. In this case, we will compare to see if height in June or September is a better predictor.

pods.reg.full_sept = glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next*herb_avg|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)

pods.reg.full_june = glmer(N_pods ~ h_apical*herb_avg + (h_apical*herb_avg|site/transect)+(h_apical*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)

aictab(list(pods.reg.full_june, pods.reg.full_sept), modnames = c('June', 'September'))
summary(pods.reg.full_sept)

Random Effects

September wins, so we will continue to use this in our model instead of height in June (h_apical). Now onto selection for random effects. We can try and drop the herb_avg term in the site/transect random effect first due to small variation.

pods.reg.mdl1 = glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next+h_apical.next:herb_avg|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
summary(pods.reg.mdl1)

#next dropping interaction term in site/transect
pods.reg.mdl2 = glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
summary(pods.reg.mdl2)

#dropping h_apical.next in site/transect
pods.reg.mdl3 = glmer(N_pods ~ h_apical.next*herb_avg + (1|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
summary(pods.reg.mdl3)

#now dropping it in year
pods.reg.mdl4 = glmer(N_pods ~ h_apical.next*herb_avg + (1|site/transect)+(herb_avg:h_apical.next+herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
summary(pods.reg.mdl4)

#dropping interaction in year
pods.reg.mdl5 = glmer(N_pods ~ h_apical.next*herb_avg + (1|site/transect)+(herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
summary(pods.reg.mdl5)

#trying dropping site/transect
pods.reg.mdl6 = glmer(N_pods ~ h_apical.next*herb_avg + (herb_avg:h_apical.next+herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
pods.reg.mdl7 = glmer(N_pods ~ h_apical.next*herb_avg + (herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)

aictab(list(pods.reg.full_sept, pods.reg.mdl1, pods.reg.mdl2, pods.reg.mdl3, pods.reg.mdl4, pods.reg.mdl5, pods.reg.mdl6, pods.reg.mdl7), modnames = c('pods.reg.full_sept', 'pods.reg.mdl1', 'pods.reg.mdl2', 'pods.reg.mdl3', 'pods.reg.mdl4', 'pods.reg.mdl5', 'pods.reg.mdl6', 'pods.reg.mdl7'))

Fixed Effects

Based on AICc, Model 2 is the best for random effects, continuing selection for fixed effects with this. Removing interaction term first, due to highest p-value.

pods.reg.mdl2a = glmer(N_pods ~ h_apical.next+herb_avg + (h_apical.next|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
summary(pods.reg.mdl2a)

#next goes the intercept
pods.reg.mdl2b = glmer(N_pods ~ h_apical.next+herb_avg -1 + (h_apical.next|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
summary(pods.reg.mdl2b)

#finally herb_avg
pods.reg.mdl2c = glmer(N_pods ~ h_apical.next -1 + (h_apical.next|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)
summary(pods.reg.mdl2c)

aictab(list(pods.reg.mdl2, pods.reg.mdl2a, pods.reg.mdl2b, pods.reg.mdl2c), modnames = c('pods.reg.mdl2', 'pods.reg.mdl2a', 'pods.reg.mdl2b', 'pods.reg.mdl2c'))

# Models b and c are indistinguishable, so we will go with 2b.

We'll go with 2b after looking at AICc.

pods.reg.final = glmer(N_pods ~ h_apical.next + herb_avg -1 + (h_apical.next|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)

Model Validation and Diagnostics

par(mfrow=c(2,2))
plot(fitted(pods.reg.final), residuals(pods.reg.final, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(pods.reg.final, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(pods.reg.final, "deviance"), xlab="Herbivory average", ylab="Deviance residuals")
par(mfrow=c(1,1))

DHARMa Diagnostics

pods.simulationOutput = DHARMa::simulateResiduals(fittedModel = pods.reg.final, n=1000)

DHARMa::plotSimulatedResiduals(pods.simulationOutput)
DHARMa::plotResiduals(metadata_sc$h_apical, pods.simulationOutput$scaledResiduals)
DHARMa::plotResiduals(metadata_sc$herb_avg, pods.simulationOutput$scaledResiduals)

NAs are being produced in DHARMa's simulatedResponse because of the way it resimulates random effects. After consulting DHARMa's issues page on GitHub here, there are two ways around this issue:

Conditional Random Effects

1) Try not resimulating the random effects (i.e. calculate residuals conditional on the fitted random effects).

pods.simulationOutput = DHARMa::simulateResiduals(fittedModel = pods.reg.final, n=1000, use.u = T)
DHARMa::plotSimulatedResiduals(pods.simulationOutput)
DHARMa::plotResiduals(metadata_sc$h_apical, pods.simulationOutput$scaledResiduals)
DHARMa::plotResiduals(metadata_sc$herb_avg, pods.simulationOutput$scaledResiduals)

This looks like a pretty good fit! For the sake of completeness, we will also try... 2) changing the link function. In our case from log to sqrt.

Sqrt Link Function

Now onto selection for random effects. We can try and drop the interaction term in the site/transect random effect first due to small variation.

pods.reg.sq1 = glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next+herb_avg|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)
summary(pods.reg.mdl1)

pods.reg.sq2 = glmer(N_pods ~ h_apical.next*herb_avg + (herb_avg|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)
summary(pods.reg.mdl2)

pods.reg.sq3 = glmer(N_pods ~ h_apical.next*herb_avg + (1|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)
summary(pods.reg.mdl3)

pods.reg.sq4 = glmer(N_pods ~ h_apical.next*herb_avg + (1|site/transect)+(h_apical.next+herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)
summary(pods.reg.mdl4)

pods.reg.sq5 = glmer(N_pods ~ h_apical.next*herb_avg + (1|site/transect)+(h_apical.next|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)
summary(pods.reg.mdl5)

pods.reg.sq6 = glmer(N_pods ~ h_apical.next*herb_avg + (1|site/transect)+(1|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)
summary(pods.reg.mdl6)

#winner from older model selection for comparison
pods.reg.sq7 = glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next+herb_avg|site/transect)+(h_apical.next+herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link = 'sqrt'), control=glmerCtrl)

aictab(list(pods.reg.full_sept, pods.reg.sq1, pods.reg.sq2, pods.reg.sq3, pods.reg.sq4, pods.reg.sq5, pods.reg.sq6, pods.reg.sq7), modnames = c('pods.reg.full_sept', 'pods.reg.sq1', 'pods.reg.sq2', 'pods.reg.sq3', 'pods.reg.sq4', 'pods.reg.sq5', 'pods.reg.sq6', 'pods.reg.sq7'))

# We'll move on to fixed effects with Model 1.

Fixed Effects.

pods.reg.sq1a = glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next+herb_avg|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)
summary(pods.reg.sq1a)

# interaction gets dropped first since highest p-value
pods.reg.sq1b = glmer(N_pods ~ h_apical.next + herb_avg + (h_apical.next+herb_avg|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)
summary(pods.reg.sq1b) #all terms are now stat. sig.

#and also compare a previous winner
pods.reg.sq1c = glmer(N_pods ~ h_apical.next + herb_avg - 1 + (h_apical.next+herb_avg|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)

aictab(list(pods.reg.sq1a, pods.reg.sq1b, pods.reg.sq1c), modnames = c('pods.reg.sq1a', 'pods.reg.sq1b', 'pods.reg.sq1c'))

# Going with Model 1b

Setting final model.

pods.reg.final.sq <- glmer(N_pods ~ h_apical.next + herb_avg + (h_apical.next+herb_avg|site/transect)+(h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='sqrt'), control=glmerCtrl)

DHARMa Diagnostics for model validation.

pods.simulationOutput.sq = DHARMa::simulateResiduals(fittedModel = pods.reg.final.sq, n=1000)

DHARMa::plotSimulatedResiduals(pods.simulationOutput.sq)
DHARMa::plotResiduals(metadata_sc$h_apical, pods.simulationOutput.sq$scaledResiduals)
DHARMa::plotResiduals(metadata_sc$herb_avg, pods.simulationOutput.sq$scaledResiduals)

These plots are not as strong as those for conditional random effects. Therefore we will keep a log link function and the winning model as our final model.

# Create model class
scaled <- list(h_apical = h_apical_sc, herb_avg = herb_avg_sc, h_apical.next = h_apical.next_sc)
pods.fit <- mwMod(list(mdl = pods.reg.final, vars = c("h_apical.next", "herb_avg"), scaled = scaled))
checkPars(pods.fit) # Check parameters
summary(pods.reg.final)

Check for overdispersion

From GLMM Wiki, we grabbed and use the overdisp_fun function.

overdisp_fun <- function(model) {
  ## number of variance parameters in 
  ##   an n-by-n variance-covariance matrix
  vpars <- function(m) {
    nrow(m)*(nrow(m)+1)/2
  }
  model.df <- sum(sapply(VarCorr(model),vpars))+length(fixef(model))
  rdf <- nrow(model.frame(model))-model.df
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

overdisp_fun(pods.reg.final)

References

Appendix

glmmLasso

The next bit is stuff from the function glmmLasso, which is another package that does variable selection for GLMMs. Keeping it for reference (took me a bit to figure out how to use the damn thing).

gL1 <- glmmLasso(fec.flower ~ h_apical + herb_avg + h_apical:herb_avg, rnd = list(site=~1, year=~1+h_apical), lambda=10, data=metadata_sc, family=binomial(link=logit))
gL2 <- glmmLasso(fec.flower ~ h_apical + herb_avg + h_apical:herb_avg, rnd = list(site=~1+herb_avg, year=~1+h_apical), lambda=10, data=metadata_sc, family=binomial(link=logit))

To use it properly, lambda needs to be optimized. Type demo("glmmLasso-soccer") in R to see how this is done (draft code is in the file glmmFit.R).

Growth (individual sites)

Let's look at BLD1 first:

metaBLD1 <- metadata %>% filter(site == "BLD1",
                                !is.na(h_apical),
                                !is.na(h_apical.next),
                                !is.na(herb_avg),
                                fec.flower == 1,
                                surv == 1)
xyplot(h_apical.next ~ h_apical, data=metaBLD1, group = year, auto.key = TRUE)

growth.reg.null=lmer(h_apical.next~h_apical*herb_avg+(1|year),data=metaBLD1, REML=T)

growth.reg.height=lmer(h_apical.next~h_apical*herb_avg+(h_apical|year),data=metaBLD1, REML=T)

growth.reg.herb=lmer(h_apical.next~h_apical*herb_avg+(herb_avg|year),data=metaBLD1, REML=T)

AIC(growth.reg.null, growth.reg.height, growth.reg.herb)

growth.mdl1.BLD1 <- lm(h_apical.next ~ h_apical*herb_avg, data=metaBLD1)
growth.mdlnoi.BLD1 <- lm(h_apical.next ~ h_apical+herb_avg, data=metaBLD1)
growth.mdlnohe.BLD1 <- lm(h_apical.next ~ h_apical+h_apical:herb_avg, data=metaBLD1)
growth.mdl0.BLD1 <- lm(h_apical.next ~ h_apical, data=metaBLD1)

AIC(growth.mdl1.BLD1, growth.mdlnoi.BLD1, growth.mdlnohe.BLD1, growth.mdl0.BLD1)

par(mfrow=c(2,2))
plot(growth.mdlnohe.BLD1)
par(mfrow=c(1,1))

Now BLD2:

metaBLD2 <- metadata %>% filter(site == "BLD2",
                                fec.flower == 1,
                                !is.na(h_apical),
                                !is.na(h_apical.next),
                                !is.na(herb_avg))
growth.mdl1.BLD2 <- lm(h_apical.next ~ h_apical*herb_avg, data=metaBLD2)
growth.mdlnoi.BLD2 <- lm(h_apical.next ~ h_apical+herb_avg, data=metaBLD2)
growth.mdlnohe.BLD2 <- lm(h_apical.next ~ h_apical+h_apical:herb_avg, data=metaBLD2)
growth.mdl0.BLD2 <- lm(h_apical.next ~ h_apical, data=metaBLD2)

AIC(growth.mdl1.BLD2, growth.mdlnoi.BLD2, growth.mdlnohe.BLD2, growth.mdl0.BLD2)

par(mfrow=c(2,2))
plot(growth.mdlnohe.BLD2)
par(mfrow=c(1,1))

Now PWR:

metaPWR <- metadata %>% filter(site == "PWR",
                                fec.flower == 1,
                                !is.na(h_apical),
                                !is.na(h_apical.next),
                                !is.na(herb_avg))
growth.mdl1.PWR <- lm(h_apical.next ~ h_apical*herb_avg, data=metaPWR)
growth.mdlnoi.PWR <- lm(h_apical.next ~ h_apical+herb_avg, data=metaPWR)
growth.mdlnohe.PWR <- lm(h_apical.next ~ h_apical+h_apical:herb_avg, data=metaPWR)
growth.mdl0.PWR <- lm(h_apical.next ~ h_apical, data=metaPWR)

AIC(growth.mdl1.PWR, growth.mdlnoi.PWR, growth.mdlnohe.PWR, growth.mdl0.PWR)

par(mfrow=c(2,2))
plot(growth.mdlnohe.PWR)
par(mfrow=c(1,1))

Now SKY:

metaSKY <- metadata %>% filter(site == "SKY",
                                fec.flower == 1,
                                !is.na(h_apical),
                                !is.na(h_apical.next),
                                !is.na(herb_avg))
growth.mdl1.SKY <- lm(h_apical.next ~ h_apical*herb_avg, data=metaSKY)
growth.mdlnoi.SKY <- lm(h_apical.next ~ h_apical+herb_avg, data=metaSKY)
growth.mdlnohe.SKY <- lm(h_apical.next ~ h_apical+h_apical:herb_avg, data=metaSKY)
growth.mdl0.SKY <- lm(h_apical.next ~ h_apical, data=metaSKY)

AIC(growth.mdl1.SKY, growth.mdlnoi.SKY, growth.mdlnohe.SKY, growth.mdl0.SKY)

par(mfrow=c(2,2))
plot(growth.mdl1.SKY)
par(mfrow=c(1,1))

Now YTB:

metaYTB <- metadata %>% filter(site == "YTB",
                               fec.flower == 1,
                               surv == 1,
                               h_apical > 50,
                               h_apical.next > 50,
                               !is.na(h_apical),
                               !is.na(h_apical.next),
                               !is.na(herb_avg))

xyplot(h_apical.next ~ h_apical, data=metaYTB, group = year, auto.key = TRUE)

growth.reg.null=lmer(h_apical.next~h_apical*herb_avg+(1|year),data=metaYTB, REML=T)

growth.reg.height=lmer(h_apical.next~h_apical*herb_avg+(h_apical|year),data=metaYTB, REML=T)

AIC(growth.reg.null, growth.reg.height)

growth.mdl1.YTB <- lm(h_apical.next ~ h_apical*herb_avg, data=metaYTB)
growth.mdlnoi.YTB <- lm(h_apical.next ~ h_apical+herb_avg, data=metaYTB)
growth.mdlnohe.YTB <- lm(h_apical.next ~ h_apical+h_apical:herb_avg, data=metaYTB)
growth.mdl0.YTB <- lm(h_apical.next ~ h_apical, data=metaYTB)

AIC(growth.mdl1.YTB, growth.mdlnoi.YTB, growth.mdlnohe.YTB, growth.mdl0.YTB)

par(mfrow=c(2,2))
plot(growth.mdl1.YTB)
par(mfrow=c(1,1))