mdlama/milkweed

Package index
Search the mdlama/milkweed package
Vignettes
Functions
89
Source code
4
Man pages
49
Home
/
GitHub
/
mdlama/milkweed
/

In mdlama/milkweed:

Initialization

library(milkweed)
library(dplyr)
library(tidyr)
library(ggplot2)
library(lme4)
library(AICcmodavg)
library(vegan)
library(flextable)
library(webshot)

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))

mwCache <- getOption("milkweed.cache")

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 = ~ 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)

Random Effects - Starting with Site

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

compute <- params$compute
savefile <- "aictabs_F_RE_Site.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 3 terms
  flower.re.site_haThe.year_haThe <- 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())

  # 2 terms
  flower.re.site_hePhaIhe.year_haThe <- glmer(fec.flower ~ h_apical*herb_avg + (herb_avg + h_apical:herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

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

  flower.re.site_haPhe.year_haThe <- 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())

  # 1 term
  flower.re.site_he.year_haThe <- glmer(fec.flower ~ h_apical*herb_avg + (herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

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

  flower.re.site_haIhe.year_haThe <- 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())

  # 0 terms
  flower.re.site_1.year_haThe <- glmer(fec.flower ~ h_apical*herb_avg + (1|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  save(flower.re.site_haThe.year_haThe,
       flower.re.site_hePhaIhe.year_haThe,
       flower.re.site_haPhaIhe.year_haThe,
       flower.re.site_haPhe.year_haThe,
       flower.re.site_he.year_haThe,
       flower.re.site_ha.year_haThe,
       flower.re.site_haIhe.year_haThe,
       flower.re.site_1.year_haThe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

aictab(list(flower.re.site_haThe.year_haThe,
            flower.re.site_hePhaIhe.year_haThe,
            flower.re.site_haPhaIhe.year_haThe,
            flower.re.site_haPhe.year_haThe,
            flower.re.site_he.year_haThe,
            flower.re.site_ha.year_haThe,
            flower.re.site_haIhe.year_haThe,
            flower.re.site_1.year_haThe),
       modnames = c("Site: height*herb", 
                    "Site: herb + height:herb", 
                    "Site: height + height:herb", 
                    "Site: height + herb",
                    "Site: herb",
                    "Site: height",
                    "Site: height:herb",
                    "Site: 1"))

Height wins! Let's move to year:

compute <- params$compute
savefile <- "aictabs_F_RE_Year.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 2 terms
  flower.re.site_ha.year_haPhe <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

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

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

  # 1 term
  flower.re.site_ha.year_he <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical|site/transect) + (herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

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

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

  # 0 terms
  flower.re.site_ha.year_1 <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical|site/transect) + (1|year), data=metadata_sc, nAGQ=0, family=binomial())

  save(flower.re.site_ha.year_haThe,
       flower.re.site_ha.year_hePhaIhe,
       flower.re.site_ha.year_haPhaIhe,
       flower.re.site_ha.year_haPhe,
       flower.re.site_ha.year_he,
       flower.re.site_ha.year_ha,
       flower.re.site_ha.year_haIhe,
       flower.re.site_ha.year_1, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

aictab(list(flower.re.site_ha.year_haThe,
            flower.re.site_ha.year_hePhaIhe,
            flower.re.site_ha.year_haPhaIhe,
            flower.re.site_ha.year_haPhe,
            flower.re.site_ha.year_he,
            flower.re.site_ha.year_ha,
            flower.re.site_ha.year_haIhe,
            flower.re.site_ha.year_1),
       modnames = c("Year: height*herb", 
                    "Year: herb + height:herb", 
                    "Year: height + height:herb", 
                    "Year: height + herb",
                    "Year: herb",
                    "Year: height",
                    "Year: height:herb",
                    "Year: 1"))

compute <- params$compute
savefile <- "aictabs_F_Sing_1.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  flower.fe.haThe.re.site_ha.year_haThe <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  save(flower.fe.haThe.re.site_ha.year_haThe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}
rePCA(flower.fe.haThe.re.site_ha.year_haThe)

We have a singularity issue. Looks like issue is in the year random effect. Let's drop the interaction term.

compute <- params$compute
savefile <- "aictabs_F_Sing_2.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  flower.fe.haThe.re.site_ha.year_haPhe <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  save(flower.fe.haThe.re.site_ha.year_haPhe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}
rePCA(flower.fe.haThe.re.site_ha.year_haPhe)

Still have a problem, but now in site. Let's drop height:

compute <- params$compute
savefile <- "aictabs_f_Sing_3.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  flower.fe.haThe.re.site_1.year_haPhe <- glmer(fec.flower ~ h_apical*herb_avg + (1|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  save(flower.fe.haThe.re.site_1.year_haPhe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}
rePCA(flower.fe.haThe.re.site_1.year_haPhe)

Random Effects - Starting with Year

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

compute <- params$compute
savefile <- "aictabs_F_RE_Year_First.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 3 terms
  flower.re.year_haThe.site_haThe <- 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())

  # 2 terms
  flower.re.year_hePhaIhe.site_haThe <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical*herb_avg|site/transect) + (herb_avg + h_apical:herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

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

  flower.re.year_haPhe.site_haThe <- 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())

  # 1 term
  flower.re.year_he.site_haThe <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical*herb_avg|site/transect) + (herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

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

  flower.re.year_haIhe.site_haThe <- 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())

  # 0 terms
  flower.re.year_1.site_haThe <- glmer(fec.flower ~ h_apical*herb_avg + (h_apical*herb_avg|site/transect) + (1|year), data=metadata_sc, nAGQ=0, family=binomial())

  save(flower.re.year_haThe.site_haThe,
       flower.re.year_hePhaIhe.site_haThe,
       flower.re.year_haPhaIhe.site_haThe,
       flower.re.year_haPhe.site_haThe,
       flower.re.year_he.site_haThe,
       flower.re.year_ha.site_haThe,
       flower.re.year_haIhe.site_haThe,
       flower.re.year_1.site_haThe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

aictab(list(flower.re.year_haThe.site_haThe,
            flower.re.year_hePhaIhe.site_haThe,
            flower.re.year_haPhaIhe.site_haThe,
            flower.re.year_haPhe.site_haThe,
            flower.re.year_he.site_haThe,
            flower.re.year_ha.site_haThe,
            flower.re.year_haIhe.site_haThe,
            flower.re.year_1.site_haThe),
       modnames = c("Year: height*herb", 
                    "Year: herb + height:herb", 
                    "Year: height + height:herb", 
                    "Year: height + herb",
                    "Year: herb",
                    "Year: height",
                    "Year: height:herb",
                    "Year: 1"))

Fixed Effects

compute <- params$compute
savefile <- "aictabs_F_FE.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 2 terms
  flower.fe.haPhe.re.site_1.year_haPhe <- glmer(fec.flower ~ h_apical+herb_avg + (1|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  flower.fe.haPhaIhe.re.site_1.year_haPhe <- glmer(fec.flower ~ h_apical+h_apical:herb_avg + (1|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  flower.fe.hePhaIhe.re.site_1.year_haPhe <- glmer(fec.flower ~ herb_avg+h_apical:herb_avg + (1|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  # 1 term
  flower.fe.ha.re.site_1.year_haPhe <- glmer(fec.flower ~ h_apical + (1|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  flower.fe.he.re.site_1.year_haPhe <- glmer(fec.flower ~ herb_avg + (1|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  flower.fe.haIhe.re.site_1.year_haPhe <- glmer(fec.flower ~ h_apical:herb_avg + (1|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  # 0 terms
  flower.fe.1.re.site_1.year_haPhe <- glmer(fec.flower ~ 1 + (1|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  save(flower.fe.haThe.re.site_1.year_haPhe,
       flower.fe.hePhaIhe.re.site_1.year_haPhe,
       flower.fe.haPhaIhe.re.site_1.year_haPhe,
       flower.fe.haPhe.re.site_1.year_haPhe,
       flower.fe.he.re.site_1.year_haPhe,
       flower.fe.ha.re.site_1.year_haPhe,
       flower.fe.haIhe.re.site_1.year_haPhe,
       flower.fe.1.re.site_1.year_haPhe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

aictab(list(flower.fe.haThe.re.site_1.year_haPhe,
            flower.fe.hePhaIhe.re.site_1.year_haPhe,
            flower.fe.haPhaIhe.re.site_1.year_haPhe,
            flower.fe.haPhe.re.site_1.year_haPhe,
            flower.fe.he.re.site_1.year_haPhe,
            flower.fe.ha.re.site_1.year_haPhe,
            flower.fe.haIhe.re.site_1.year_haPhe,
            flower.fe.1.re.site_1.year_haPhe),
       modnames = c("Fixed: height*herb", 
                    "Fixed: herb + height:herb", 
                    "Fixed: height + height:herb", 
                    "Fixed: height + herb",
                    "Fixed: herb",
                    "Fixed: height",
                    "Fixed: height:herb",
                    "Fixed: 1"))

Winner:

summary(flower.fe.haThe.re.site_1.year_haPhe)

Model 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.fe.haThe.re.site_1.year_haPhe), residuals(flower.fe.haThe.re.site_1.year_haPhe, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(flower.fe.haThe.re.site_1.year_haPhe, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(flower.fe.haThe.re.site_1.year_haPhe, "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.fe.haThe.re.site_1.year_haPhe, 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.

Create AIC Tables

Random Effect: Site

Let's create the full table for publication:

f.re.site.AIC <- aictab(list(flower.re.site_haThe.year_haThe,
            flower.re.site_hePhaIhe.year_haThe,
            flower.re.site_haPhaIhe.year_haThe,
            flower.re.site_haPhe.year_haThe,
            flower.re.site_he.year_haThe,
            flower.re.site_ha.year_haThe,
            flower.re.site_haIhe.year_haThe,
            flower.re.site_1.year_haThe),
       modnames = c("Site: height*herb", 
                    "Site: herb + height:herb", 
                    "Site: height + height:herb", 
                    "Site: height + herb",
                    "Site: herb",
                    "Site: height",
                    "Site: height:herb",
                    "Site: 1"))
f.re.site.AIC

Let's format this puppy:

ft <- flextable(f.re.site.AIC) %>%
  add_header(Modnames = "Flowering: Site Random Effect") %>%
  fontsize(part = "header", size = 12) %>%
  set_header_labels(Modnames = "Model") %>% 
  bg(i = 5, bg = "#dedede") %>% 
  width(j = ~ Modnames, width = 2.6) %>%
  set_formatter(
    K = function(x) sprintf("%d", x)
  )
ft

Now save as image:

# create an Rmd file ----
rmd_name <- tempfile(fileext = ".Rmd")
cat("```r\nft\n```", file = rmd_name)
# render as an html file ----
html_name <- tempfile(fileext = ".html")
rmarkdown::render(rmd_name, output_format = "html_document", output_file = html_name)
# get a png from the html file with webshot ----
webshot(html_name, zoom = 2, file = "floweringRE_Site_AIC.png", 
        selector = "body > div.container-fluid.main-container > div.tabwid > table")

Random Effect: Year

Let's create the full table for publication:

f.re.year.AIC <- aictab(list(flower.re.site_ha.year_haThe,
            flower.re.site_ha.year_hePhaIhe,
            flower.re.site_ha.year_haPhaIhe,
            flower.re.site_ha.year_haPhe,
            flower.re.site_ha.year_he,
            flower.re.site_ha.year_ha,
            flower.re.site_ha.year_haIhe,
            flower.re.site_ha.year_1),
       modnames = c("Year: height*herb", 
                    "Year: herb + height:herb", 
                    "Year: height + height:herb", 
                    "Year: height + herb",
                    "Year: herb",
                    "Year: height",
                    "Year: height:herb",
                    "Year: 1"))
f.re.year.AIC

Let's format this puppy:

ft <- flextable(f.re.year.AIC) %>%
  add_header(Modnames = "Flowering: Year Random Effect") %>%
  fontsize(part = "header", size = 12) %>%
  set_header_labels(Modnames = "Model") %>% 
  bg(i = 3, bg = "#dedede") %>% 
  width(j = ~ Modnames, width = 2.6) %>%
  set_formatter(
    K = function(x) sprintf("%d", x)
  )
ft

Now save as image:

# create an Rmd file ----
rmd_name <- tempfile(fileext = ".Rmd")
cat("```r\nft\n```", file = rmd_name)
# render as an html file ----
html_name <- tempfile(fileext = ".html")
rmarkdown::render(rmd_name, output_format = "html_document", output_file = html_name)
# get a png from the html file with webshot ----
webshot(html_name, zoom = 2, file = "floweringRE_Year_AIC.png", 
        selector = "body > div.container-fluid.main-container > div.tabwid > table")

Dealing with singularity

f.singularity.AIC <- aictab(list(flower.fe.haThe.re.site_ha.year_haThe,
                                 flower.fe.haThe.re.site_ha.year_haPhe,
                                 flower.fe.haThe.re.site_1.year_haPhe),
                            modnames = c("(S) Site: height, Year: h_apical*height",
                                         "(S) Site: height, Year: h_apical+height",
                                         "(NS) Site: 1, Year: h_apical+height"))
f.singularity.AIC

Let's format this puppy:

ft <- flextable(f.singularity.AIC) %>%
  add_header(Modnames = "Flowering: Singularity Adjustments") %>%
  fontsize(part = "header", size = 12) %>%
  set_header_labels(Modnames = "Model") %>% 
  bg(i = 3, bg = "#dedede") %>% 
  width(j = ~ Modnames, width = 2.6) %>%
  set_formatter(
    K = function(x) sprintf("%d", x)
  )
ft

Now save as image:

# create an Rmd file ----
rmd_name <- tempfile(fileext = ".Rmd")
cat("```r\nft\n```", file = rmd_name)
# render as an html file ----
html_name <- tempfile(fileext = ".html")
rmarkdown::render(rmd_name, output_format = "html_document", output_file = html_name)
# get a png from the html file with webshot ----
# webshot(html_name, zoom = 2, file = "floweringRE_Singularity_AIC.png", 
#         selector = "body > div.container-fluid.main-container > div.tabwid > table")

Fixed Effects

Let's create the full table for publication:

f.fe.AIC <- aictab(list(flower.fe.haThe.re.site_1.year_haPhe,
            flower.fe.hePhaIhe.re.site_1.year_haPhe,
            flower.fe.haPhaIhe.re.site_1.year_haPhe,
            flower.fe.haPhe.re.site_1.year_haPhe,
            flower.fe.he.re.site_1.year_haPhe,
            flower.fe.ha.re.site_1.year_haPhe,
            flower.fe.haIhe.re.site_1.year_haPhe,
            flower.fe.1.re.site_1.year_haPhe),
       modnames = c("Fixed: height*herb", 
                    "Fixed: herb + height:herb", 
                    "Fixed: height + height:herb", 
                    "Fixed: height + herb",
                    "Fixed: herb",
                    "Fixed: height",
                    "Fixed: height:herb",
                    "Fixed: 1"))
f.fe.AIC

Let's format this puppy:

ft <- flextable(f.fe.AIC) %>%
  add_header(Modnames = "Flowering: Fixed Effects") %>%
  fontsize(part = "header", size = 12) %>%
  set_header_labels(Modnames = "Model") %>% 
  bg(i = 2, bg = "#dedede") %>% 
  width(j = ~ Modnames, width = 2.6) %>%
  set_formatter(
    K = function(x) sprintf("%d", x)
  )
ft

Now save as image:

# create an Rmd file ----
rmd_name <- tempfile(fileext = ".Rmd")
cat("```r\nft\n```", file = rmd_name)
# render as an html file ----
html_name <- tempfile(fileext = ".html")
rmarkdown::render(rmd_name, output_format = "html_document", output_file = html_name)
# get a png from the html file with webshot ----
# webshot(html_name, zoom = 2, file = "floweringFE_AIC.png", 
#         selector = "body > div.container-fluid.main-container > div.tabwid > table")

Visualizations

Create model class

scaled <- list(h_apical = h_apical_sc, herb_avg = herb_avg_sc)
flower.fit <- mwMod(list(mdl = flower.fe.haThe.re.site_1.year_haPhe, vars = c("h_apical", "herb_avg"), scaled = scaled))
checkPars(flower.fit) # Check parameters

flower.fit$pars$scaled

Unscaled parameters across site:

flower.fit$pars$unscaled

Let's plot full data.

# 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 <- 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 = ~ 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)

Random Effects

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

compute <- params$compute
savefile <- "aictabs_S_RE_Site.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 3 terms
  surv.re.site_haThe.year_haThe <- glmer(surv ~ h_apical*herb_avg + (h_apical*herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  # 2 terms
  surv.re.site_hePhaIhe.year_haThe <- glmer(surv ~ h_apical*herb_avg + (herb_avg + h_apical:herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  surv.re.site_haPhaIhe.year_haThe <- glmer(surv ~ h_apical*herb_avg + (h_apical + h_apical:herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  surv.re.site_haPhe.year_haThe <- glmer(surv ~ h_apical*herb_avg + (h_apical + herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  # 1 term
  surv.re.site_he.year_haThe <- glmer(surv ~ h_apical*herb_avg + (herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  surv.re.site_ha.year_haThe <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  surv.re.site_haIhe.year_haThe <- glmer(surv ~ h_apical*herb_avg + (h_apical:herb_avg|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  # 0 terms
  surv.re.site_1.year_haThe <- glmer(surv ~ h_apical*herb_avg + (1|site/transect) + (h_apical*herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  save(surv.re.site_haThe.year_haThe,
       surv.re.site_hePhaIhe.year_haThe,
       surv.re.site_haPhaIhe.year_haThe,
       surv.re.site_haPhe.year_haThe,
       surv.re.site_he.year_haThe,
       surv.re.site_ha.year_haThe,
       surv.re.site_haIhe.year_haThe,
       surv.re.site_1.year_haThe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

s.re.site.AIC <- aictab(list(surv.re.site_haThe.year_haThe,
            surv.re.site_hePhaIhe.year_haThe,
            surv.re.site_haPhaIhe.year_haThe,
            surv.re.site_haPhe.year_haThe,
            surv.re.site_he.year_haThe,
            surv.re.site_ha.year_haThe,
            surv.re.site_haIhe.year_haThe,
            surv.re.site_1.year_haThe),
       modnames = c("Site: height*herb", 
                    "Site: herb + height:herb", 
                    "Site: height + height:herb", 
                    "Site: height + herb",
                    "Site: herb",
                    "Site: height",
                    "Site: height:herb",
                    "Site: 1"))
s.re.site.AIC

Height wins. Let's move to year:

compute <- params$compute
savefile <- "aictabs_S_RE_Year.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 2 terms
  surv.re.site_ha.year_haPhe <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  surv.re.site_ha.year_haPhaIhe <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical+h_apical:herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  surv.re.site_ha.year_hePhaIhe <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (herb_avg+h_apical:herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  # 1 term
  surv.re.site_ha.year_he <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  surv.re.site_ha.year_ha <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=0, family=binomial())

  surv.re.site_ha.year_haIhe <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical:herb_avg|year), data=metadata_sc, nAGQ=0, family=binomial())

  # 0 terms
  surv.re.site_ha.year_1 <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (1|year), data=metadata_sc, nAGQ=0, family=binomial())

  save(surv.re.site_ha.year_haThe,
       surv.re.site_ha.year_hePhaIhe,
       surv.re.site_ha.year_haPhaIhe,
       surv.re.site_ha.year_haPhe,
       surv.re.site_ha.year_he,
       surv.re.site_ha.year_ha,
       surv.re.site_ha.year_haIhe,
       surv.re.site_ha.year_1, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

s.re.year.AIC <- aictab(list(surv.re.site_ha.year_haThe,
            surv.re.site_ha.year_hePhaIhe,
            surv.re.site_ha.year_haPhaIhe,
            surv.re.site_ha.year_haPhe,
            surv.re.site_ha.year_he,
            surv.re.site_ha.year_ha,
            surv.re.site_ha.year_haIhe,
            surv.re.site_ha.year_1),
       modnames = c("Year: height*herb", 
                    "Year: herb + height:herb", 
                    "Year: height + height:herb", 
                    "Year: height + herb",
                    "Year: herb",
                    "Year: height",
                    "Year: height:herb",
                    "Year: 1"))

s.re.year.AIC

compute <- params$compute
savefile <- "aictabs_S_Sing_1.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  surv.fe.haThe.re.site_ha.year_haPhe <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical+herb_avg|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  save(surv.fe.haThe.re.site_ha.year_haPhe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}
rePCA(surv.fe.haThe.re.site_ha.year_haPhe)

We have a singularity issue. Looks like issue is in the year random effect. Let's drop herbivory.

compute <- params$compute
savefile <- "aictabs_S_Sing_2.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  surv.fe.haThe.re.site_ha.year_ha <- glmer(surv ~ h_apical*herb_avg + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  save(surv.fe.haThe.re.site_ha.year_ha, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

rePCA(surv.fe.haThe.re.site_ha.year_ha)

Fixed Effects

compute <- params$compute
savefile <- "aictabs_S_FE.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 2 terms
  surv.fe.haPhe.re.site_ha.year_ha <- glmer(surv ~ h_apical+herb_avg + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  surv.fe.haPhaIhe.re.site_ha.year_ha <- glmer(surv ~ h_apical+h_apical:herb_avg + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  surv.fe.hePhaIhe.re.site_ha.year_ha <- glmer(surv ~ herb_avg+h_apical:herb_avg + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  # 1 term
  surv.fe.ha.re.site_ha.year_ha <- glmer(surv ~ h_apical + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  surv.fe.he.re.site_ha.year_ha <- glmer(surv ~ herb_avg + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  surv.fe.haIhe.re.site_ha.year_ha <- glmer(surv ~ h_apical:herb_avg + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  # 0 terms
  surv.fe.1.re.site_ha.year_ha <- glmer(surv ~ 1 + (h_apical|site/transect) + (h_apical|year), data=metadata_sc, nAGQ=1, family=binomial(), control=glmerCtrl)

  save(surv.fe.haThe.re.site_ha.year_ha,
       surv.fe.hePhaIhe.re.site_ha.year_ha,
       surv.fe.haPhaIhe.re.site_ha.year_ha,
       surv.fe.haPhe.re.site_ha.year_ha,
       surv.fe.he.re.site_ha.year_ha,
       surv.fe.ha.re.site_ha.year_ha,
       surv.fe.haIhe.re.site_ha.year_ha,
       surv.fe.1.re.site_ha.year_ha, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

s.fe.AIC <- aictab(list(surv.fe.haThe.re.site_ha.year_ha,
            surv.fe.hePhaIhe.re.site_ha.year_ha,
            surv.fe.haPhaIhe.re.site_ha.year_ha,
            surv.fe.haPhe.re.site_ha.year_ha,
            surv.fe.he.re.site_ha.year_ha,
            surv.fe.ha.re.site_ha.year_ha,
            surv.fe.haIhe.re.site_ha.year_ha,
            surv.fe.1.re.site_ha.year_ha),
       modnames = c("Fixed: height*herb", 
                    "Fixed: herb + height:herb", 
                    "Fixed: height + height:herb", 
                    "Fixed: height + herb",
                    "Fixed: herb",
                    "Fixed: height",
                    "Fixed: height:herb",
                    "Fixed: 1"))

s.fe.AIC

Winner:

summary(surv.fe.haPhe.re.site_ha.year_ha)

Model Validation and Diagnostics

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

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

DHARMa Diagnostics. Same process as above. 

surv.simulationOutput = DHARMa::simulateResiduals(fittedModel = surv.fe.haPhe.re.site_ha.year_ha, 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 do not pass the Uniformity test, although the p-value is notably large.

DHARMa::testUniformity(surv.simulationOutput)

Create AIC Tables

Random Effect: Site

Let's format this puppy:

ft <- flextable(s.re.site.AIC) %>%
  add_header(Modnames = "Survival: Site Random Effect") %>%
  fontsize(part = "header", size = 12) %>%
  set_header_labels(Modnames = "Model") %>% 
  bg(i = 1, bg = "#dedede") %>% 
  width(j = ~ Modnames, width = 2.6) %>%
  set_formatter(
    K = function(x) sprintf("%d", x)
  )
ft

Now save as image:

# create an Rmd file ----
rmd_name <- tempfile(fileext = ".Rmd")
cat("```r\nft\n```", file = rmd_name)
# render as an html file ----
html_name <- tempfile(fileext = ".html")
rmarkdown::render(rmd_name, output_format = "html_document", output_file = html_name)
# get a png from the html file with webshot ----
# webshot(html_name, zoom = 2, file = "survivalRE_Site_AIC.png", 
#         selector = "body > div.container-fluid.main-container > div.tabwid > table")

Random Effect: Year

Let's format this puppy:

ft <- flextable(s.re.year.AIC) %>%
  add_header(Modnames = "Survival: Year Random Effect") %>%
  fontsize(part = "header", size = 12) %>%
  set_header_labels(Modnames = "Model") %>% 
  bg(i = 2, bg = "#dedede") %>% 
  width(j = ~ Modnames, width = 2.6) %>%
  set_formatter(
    K = function(x) sprintf("%d", x)
  )
ft

Now save as image:

# create an Rmd file ----
rmd_name <- tempfile(fileext = ".Rmd")
cat("```r\nft\n```", file = rmd_name)
# render as an html file ----
html_name <- tempfile(fileext = ".html")
rmarkdown::render(rmd_name, output_format = "html_document", output_file = html_name)
# get a png from the html file with webshot ----
# webshot(html_name, zoom = 2, file = "survivalRE_Year_AIC.png", 
#         selector = "body > div.container-fluid.main-container > div.tabwid > table")

Dealing with singularity

s.singularity.AIC <- aictab(list(surv.fe.haThe.re.site_ha.year_haPhe,
                                 surv.fe.haThe.re.site_ha.year_ha),
                            modnames = c("(S) Site: h_apical, Year: h_apical+height",
                                         "(NS) Site: h_apical, Year: h_apical"))
s.singularity.AIC

Let's format this puppy:

ft <- flextable(s.singularity.AIC) %>%
  add_header(Modnames = "Survival: Singularity Adjustments") %>%
  fontsize(part = "header", size = 12) %>%
  set_header_labels(Modnames = "Model") %>% 
  bg(i = 2, bg = "#dedede") %>% 
  width(j = ~ Modnames, width = 2.6) %>%
  set_formatter(
    K = function(x) sprintf("%d", x)
  )
ft

Now save as image:

# create an Rmd file ----
rmd_name <- tempfile(fileext = ".Rmd")
cat("```r\nft\n```", file = rmd_name)
# render as an html file ----
html_name <- tempfile(fileext = ".html")
rmarkdown::render(rmd_name, output_format = "html_document", output_file = html_name)
# get a png from the html file with webshot ----
# webshot(html_name, zoom = 2, file = "survivalRE_Singularity_AIC.png", 
#         selector = "body > div.container-fluid.main-container > div.tabwid > table")

Fixed Effects

Let's format this puppy:

ft <- flextable(s.fe.AIC) %>%
  add_header(Modnames = "Survival: Fixed Effects") %>%
  fontsize(part = "header", size = 12) %>%
  set_header_labels(Modnames = "Model") %>% 
  bg(i = 5, bg = "#dedede") %>% 
  width(j = ~ Modnames, width = 2.6) %>%
  set_formatter(
    K = function(x) sprintf("%d", x)
  )
ft

Now save as image:

# create an Rmd file ----
rmd_name <- tempfile(fileext = ".Rmd")
cat("```r\nft\n```", file = rmd_name)
# render as an html file ----
html_name <- tempfile(fileext = ".html")
rmarkdown::render(rmd_name, output_format = "html_document", output_file = html_name)
# get a png from the html file with webshot ----
# webshot(html_name, zoom = 2, file = "survivalFE_AIC.png", 
#         selector = "body > div.container-fluid.main-container > div.tabwid > table")

Visualizations

Create model class

# Create model class
scaled <- list(h_apical = h_apical_sc, herb_avg = herb_avg_sc)
surv.fit <- mwMod(list(mdl = surv.fe.haPhe.re.site_ha.year_ha, vars = c("h_apical", "herb_avg"), scaled = scaled))
checkPars(surv.fit) # Check parameters

surv.fit$pars$scaled

Unscaled parameters across site:

surv.fit$pars$unscaled

(Plot later)

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)

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

Random Effects

growth.re.site_haThe.year_haThe <- 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.re.site_haThe.year_haThe)

First drop interaction term on site.

growth.re.site_haPhe.year_haThe <- 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.re.site_haPhe.year_haThe)

Now drop herbivory on site.

growth.re.site_ha.year_haThe <- 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.re.site_ha.year_haThe)

Now drop interaction term on year.

growth.re.site_ha.year_haPhe <- 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.re.site_ha.year_haPhe)

Now drop herbivory on year.

growth.re.site_ha.year_ha <- lmer(h_apical.next~h_apical*herb_avg + (h_apical|site/transect)+(h_apical|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.re.site_ha.year_ha)

growth.re.transect_ha.year_haPhe <- lmer(h_apical.next~h_apical*herb_avg + (h_apical|transect)+(h_apical+herb_avg|year),data=metadata_sc, REML=T, control=lmerCtrl)

Let's check model selection with candidates.

growth.re.AIC <- aictab(list(growth.re.site_haThe.year_haThe,
            growth.re.site_haPhe.year_haThe,
            growth.re.site_ha.year_haThe,
            growth.re.site_ha.year_haPhe,
            growth.re.site_ha.year_ha,
            growth.re.transect_ha.year_haPhe),
       modnames = c('(S) Site: herbivory*height, Year: herbivory*height', 
                    '(S) Site: herbivory+height, Year: herbivory*height', 
                    '(S) Site: height, Year: herbivory*height', 
                    '(NS) Site: height, Year: herbivory+height', 
                    '(NS) Site: height, Year: height',
                    '(NS) Transect: height, Year: height+herbivory'))
sapply(list(growth.re.site_haThe.year_haThe,
            growth.re.site_haPhe.year_haThe,
            growth.re.site_ha.year_haThe,
            growth.re.site_ha.year_haPhe,
            growth.re.site_ha.year_ha,
            growth.re.transect_ha.year_haPhe), function(x) { isSingular(x) }, simplify = TRUE)
growth.re.AIC

Fixed effects

growth.fe.haThe.site_ha.year_haPhe <- lmer(h_apical.next ~ h_apical*herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

# 2 terms
growth.fe.haPhe.site_ha.year_haPhe <- lmer(h_apical.next ~ h_apical+herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

growth.fe.haPhaIhe.site_ha.year_haPhe <- lmer(h_apical.next ~ h_apical+h_apical:herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

growth.fe.hePhaIhe.site_ha.year_haPhe <- lmer(h_apical.next ~ herb_avg+h_apical:herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

# 1 term
growth.fe.ha.site_ha.year_haPhe <- lmer(h_apical.next ~ h_apical + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

growth.fe.he.site_ha.year_haPhe <- lmer(h_apical.next ~ herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

growth.fe.haIhe.site_ha.year_haPhe <- lmer(h_apical.next ~ h_apical:herb_avg + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

# 0 terms
growth.fe.1.site_ha.year_haPhe <- lmer(h_apical.next ~ 1 + (h_apical|site/transect)+(h_apical+herb_avg|year), data=metadata_sc, REML=F, control=lmerCtrl)

growth.fe.AIC <- aictab(list(growth.fe.haThe.site_ha.year_haPhe,
                             growth.fe.haPhe.site_ha.year_haPhe,
                             growth.fe.haPhaIhe.site_ha.year_haPhe,
                             growth.fe.hePhaIhe.site_ha.year_haPhe,
                             growth.fe.ha.site_ha.year_haPhe,
                             growth.fe.he.site_ha.year_haPhe,
                             growth.fe.haIhe.site_ha.year_haPhe,
                             growth.fe.1.site_ha.year_haPhe),
                        modnames = c("Fixed: height*herb",
                                     "Fixed: height+herb",
                                     "Fixed: height+height:herb",
                                     "Fixed: herb+height:herb",
                                     "Fixed: height",
                                     "Fixed: herb",
                                     "Fixed: height:herb",
                                     "Fixed: 1"))
growth.fe.AIC

Winner:

summary(growth.fe.haThe.site_ha.year_haPhe)

Model Validation and Diagnostics

Now some diagnostics (aka model validation).

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

DHARMa Diagnostics

growth.simulationOutput = DHARMa::simulateResiduals(fittedModel = growth.fe.haThe.site_ha.year_haPhe, 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.fe.haThe.site_ha.year_haPhe.log = lmer(h_apical.next ~ h_apical*log_herb_avg + (h_apical|site/transect)+(h_apical*log_herb_avg|year), data=metadata_sc, REML=F)
summary(growth.fe.haThe.site_ha.year_haPhe.log)
growth.fe.haThe.site_ha.year_haPhe.arcsin = lmer(h_apical.next ~ h_apical*arcsin_herb_avg + (h_apical|site/transect)+(h_apical*arcsin_herb_avg|year), data=metadata_sc, REML=F)
summary(growth.fe.haThe.site_ha.year_haPhe.arcsin)

#both transformations produce a much lower AIC, with log coming out on top
aictab(list(growth.fe.haThe.site_ha.year_haPhe, growth.fe.haThe.site_ha.year_haPhe.log, growth.fe.haThe.site_ha.year_haPhe.arcsin), modnames = c('growth.fe.haThe.site_ha.year_haPhe', 'growth.fe.haThe.site_ha.year_haPhe.log', 'growth.fe.haThe.site_ha.year_haPhe.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.fe.haThe.site_ha.year_haPhe.log), residuals(growth.fe.haThe.site_ha.year_haPhe.log, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(growth.fe.haThe.site_ha.year_haPhe.log, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(growth.fe.haThe.site_ha.year_haPhe.log, "deviance"), xlab="log(Herbivory average)", ylab="Deviance residuals")
par(mfrow=c(1,1))

par(mfrow=c(2,2))
plot(fitted(growth.fe.haThe.site_ha.year_haPhe.arcsin), residuals(growth.fe.haThe.site_ha.year_haPhe.arcsin, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(growth.fe.haThe.site_ha.year_haPhe.arcsin, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(growth.fe.haThe.site_ha.year_haPhe.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.fe.haThe.site_ha.year_haPhe.log, n=1000)
growth.simulationOutput.arcsin = DHARMa::simulateResiduals(fittedModel = growth.fe.haThe.site_ha.year_haPhe.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)

DHARMa::testUniformity(growth.simulationOutput.log)
DHARMa::testUniformity(growth.simulationOutput.arcsin)

# 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.fe.haThe.site_ha.year_haPhe, vars = c("h_apical", "herb_avg", "h_apical.next"), scaled = scaled))
checkPars(growth.fit) # Check parameters

Pods

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 = ~ 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)

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.

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

compute <- params$compute
savefile <- "aictabs_P_RE_Site.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 3 terms
  pods.re.site_haThe.year_haThe <- 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=0, family=poisson(link='log'))

  # 2 terms
  pods.re.site_hePhaIhe.year_haThe <- glmer(N_pods ~ h_apical.next*herb_avg + (herb_avg + h_apical.next:herb_avg|site/transect) + (h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  pods.re.site_haPhaIhe.year_haThe <- 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=0, family=poisson(link='log'))

  pods.re.site_haPhe.year_haThe <- 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=0, family=poisson(link='log'))

  # 1 term
  pods.re.site_he.year_haThe <- glmer(N_pods ~ h_apical.next*herb_avg + (herb_avg|site/transect) + (h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  pods.re.site_ha.year_haThe <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  pods.re.site_haIhe.year_haThe <- 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=0, family=poisson(link='log'))

  # 0 terms
  pods.re.site_1.year_haThe <- glmer(N_pods ~ h_apical.next*herb_avg + (1|site/transect) + (h_apical.next*herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  save(pods.re.site_haThe.year_haThe,
       pods.re.site_hePhaIhe.year_haThe,
       pods.re.site_haPhaIhe.year_haThe,
       pods.re.site_haPhe.year_haThe,
       pods.re.site_he.year_haThe,
       pods.re.site_ha.year_haThe,
       pods.re.site_haIhe.year_haThe,
       pods.re.site_1.year_haThe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

pods.re.site.AIC <- aictab(list(pods.re.site_haThe.year_haThe,
            pods.re.site_hePhaIhe.year_haThe,
            pods.re.site_haPhaIhe.year_haThe,
            pods.re.site_haPhe.year_haThe,
            pods.re.site_he.year_haThe,
            pods.re.site_ha.year_haThe,
            pods.re.site_haIhe.year_haThe,
            pods.re.site_1.year_haThe),
       modnames = c("Site: height*herb", 
                    "Site: herb + height:herb", 
                    "Site: height + height:herb", 
                    "Site: height + herb",
                    "Site: herb",
                    "Site: height",
                    "Site: height:herb",
                    "Site: 1"))
pods.re.site.AIC

Height on site wins. Let's move on the year.

compute <- params$compute
savefile <- "aictabs_P_RE_Year.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 2 terms
  pods.re.site_ha.year_haPhe <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (h_apical.next+herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  pods.re.site_ha.year_haPhaIhe <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (h_apical.next+h_apical.next:herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  pods.re.site_ha.year_hePhaIhe <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (herb_avg+h_apical.next:herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  # 1 term
  pods.re.site_ha.year_he <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  pods.re.site_ha.year_ha <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (h_apical.next|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  pods.re.site_ha.year_haIhe <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (h_apical.next:herb_avg|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  # 0 terms
  pods.re.site_ha.year_1 <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (1|year), data=metadata_sc, nAGQ=0, family=poisson(link='log'))

  save(pods.re.site_ha.year_haThe,
       pods.re.site_ha.year_hePhaIhe,
       pods.re.site_ha.year_haPhaIhe,
       pods.re.site_ha.year_haPhe,
       pods.re.site_ha.year_he,
       pods.re.site_ha.year_ha,
       pods.re.site_ha.year_haIhe,
       pods.re.site_ha.year_1, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

p.re.year.AIC <- aictab(list(pods.re.site_haPhe.year_haThe,
            pods.re.site_ha.year_hePhaIhe,
            pods.re.site_ha.year_haPhaIhe,
            pods.re.site_ha.year_haPhe,
            pods.re.site_ha.year_he,
            pods.re.site_ha.year_ha,
            pods.re.site_ha.year_haIhe,
            pods.re.site_ha.year_1),
       modnames = c("Year: height*herb", 
                    "Year: herb + height:herb", 
                    "Year: height + height:herb", 
                    "Year: height + herb",
                    "Year: herb",
                    "Year: height",
                    "Year: height:herb",
                    "Year: 1"))

p.re.year.AIC

Herbivory + height:herb wins. Let's check singularity:

compute <- params$compute
savefile <- "aictabs_P_Sing_1.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  pods.fe.haThe.re.site_ha.year_hePhaIhe <- glmer(N_pods ~ h_apical.next*herb_avg + (h_apical.next|site/transect) + (herb_avg+h_apical.next:herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'))

  save(pods.fe.haThe.re.site_ha.year_hePhaIhe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}
rePCA(pods.fe.haThe.re.site_ha.year_hePhaIhe)

Because the model failed to converge, it looks singular. Let's try dropping the variable in site.

compute <- params$compute
savefile <- "aictabs_P_Sing_2.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  pods.fe.haThe.re.site_1.year_hePhaIhe <- 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'))

  save(pods.fe.haThe.re.site_1.year_hePhaIhe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}
rePCA(pods.fe.haThe.re.site_1.year_hePhaIhe)

This fixed the singularity in our model!

Fixed effects

compute <- params$compute
savefile <- "aictabs_P_FE.RData"

if (!file.exists(file.path(mwCache,savefile)) || (compute)) {
  # 2 terms
  pods.fe.haPhe.re.site_1.year_hePhaIhe <- 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)

  pods.fe.haPhaIhe.re.site_1.year_hePhaIhe <- glmer(N_pods ~ h_apical.next+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)

  pods.fe.hePhaIhe.re.site_1.year_hePhaIhe <- glmer(N_pods ~ herb_avg+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)

  # 1 term
  pods.fe.ha.re.site_1.year_hePhaIhe <- glmer(N_pods ~ h_apical.next + (1|site/transect) + (herb_avg+h_apical.next:herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)

  pods.fe.he.re.site_1.year_hePhaIhe <- glmer(N_pods ~ herb_avg + (1|site/transect) + (herb_avg+h_apical.next:herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)

  pods.fe.haIhe.re.site_1.year_hePhaIhe <- 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)

  # 0 terms
  pods.fe.1.re.site_1.year_hePhaIhe <- glmer(N_pods ~ 1 + (1|site/transect) + (herb_avg+h_apical.next:herb_avg|year), data=metadata_sc, nAGQ=1, family=poisson(link='log'), control=glmerCtrl)

  save(pods.fe.haThe.re.site_1.year_hePhaIhe,
       pods.fe.hePhaIhe.re.site_1.year_hePhaIhe,
       pods.fe.haPhaIhe.re.site_1.year_hePhaIhe,
       pods.fe.haPhe.re.site_1.year_hePhaIhe,
       pods.fe.he.re.site_1.year_hePhaIhe,
       pods.fe.ha.re.site_1.year_hePhaIhe,
       pods.fe.haIhe.re.site_1.year_hePhaIhe,
       pods.fe.1.re.site_1.year_hePhaIhe, 
       file=file.path(mwCache,savefile))
} else {
  load(file.path(mwCache,savefile))
}

p.fe.AIC <- aictab(list(pods.fe.haThe.re.site_1.year_hePhaIhe,
            pods.fe.hePhaIhe.re.site_1.year_hePhaIhe,
            pods.fe.haPhaIhe.re.site_1.year_hePhaIhe,
            pods.fe.haPhe.re.site_1.year_hePhaIhe,
            pods.fe.he.re.site_1.year_hePhaIhe,
            pods.fe.ha.re.site_1.year_hePhaIhe,
            pods.fe.haIhe.re.site_1.year_hePhaIhe,
            pods.fe.1.re.site_1.year_hePhaIhe),
       modnames = c("Fixed: height*herb", 
                    "Fixed: herb + height:herb", 
                    "Fixed: height + height:herb", 
                    "Fixed: height + herb",
                    "Fixed: herb",
                    "Fixed: height",
                    "Fixed: height:herb",
                    "Fixed: 1"))
p.fe.AIC

Height and herbivory is within two AIC points of the winning model, so we will use that one.

Winner:

summary(pods.fe.haPhe.re.site_1.year_hePhaIhe)

Model Validation and Diagnostics

par(mfrow=c(2,2))
plot(fitted(pods.fe.haPhe.re.site_1.year_hePhaIhe), residuals(pods.fe.haPhe.re.site_1.year_hePhaIhe, "deviance"), xlab="Fitted values", ylab="Deviance residuals")
plot(metadata_sc$h_apical, residuals(pods.fe.haThe.re.site_1.year_hePhaIhe, "deviance"), xlab="Apical height", ylab="Deviance residuals")
plot(metadata_sc$herb_avg, residuals(pods.fe.haThe.re.site_1.year_hePhaIhe, "deviance"), xlab="Herbivory average", ylab="Deviance residuals")
par(mfrow=c(1,1))

DHARMa Diagnostics

pods.simulationOutput = DHARMa::simulateResiduals(fittedModel = pods.fe.haPhe.re.site_1.year_hePhaIhe, n=1000)

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


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

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close