In karinorman/biodivTS: What the Package Does (One Line, Title Case)

library(tidyverse)
library(lme4)
library(broom.mixed)
library(lmerTest)
library(emmeans)

load(here::here("data/model_data.rda"))
amph_data <- model_data %>%
  filter(taxa == "Amphibians")

Model fitting function

get_model <- function(metric_name, formula, model_id, data, ...){
  print(c(metric_name, model_id))
  qlmer <- quietly(lmer)

  fit_data <- data %>%
    filter(metric == metric_name) %>%
    mutate(year_scaled = scale(year),
           duration_scaled = scale(duration, scale = FALSE),
           startyear_scaled = scale(startyear, scale = FALSE),
           taxa = as.factor(taxa),
           realm = as.factor(realm),
           climate = as.factor(climate))

  if(n_distinct(fit_data$rarefyID) < 2){
    return(tibble())
  }

  if(!model_id %in% c("base", "base_dur", "base_srt")){
    if(fit_data %>% pull(model_id) %>% n_distinct() >1) {

      fit <- qlmer(formula = formula, data = fit_data, ...)

      slopes <- get_fixed_slopes(metric = metric_name, model = fit$result, model_id = model_id, model_data = fit_data) 

      group_slopes <- as_tibble(slopes$emtrends) %>%
        mutate(effect = "group_slope", group = model_id) %>%
        rename(std.error = SE, estimate = year_scaled.trend, term = model_id)

      slope_constrasts <- as_tibble(slopes$contrasts) %>%
        mutate(effect = "contrast", group = model_id) %>%
        rename(std.error = SE, term = contrast, statistic = t.ratio)

      df <- tidy(fit$result) %>% 
        bind_rows(group_slopes, slope_constrasts) %>%
        mutate(metric = metric_name, model = list(fit$result), model_id = model_id,
               n_ts = n_distinct(fit_data$rarefyID),
               formula = formula,
               warning = ifelse(length(fit$warnings) == 0, NA, fit$warnings))
    }else{
      return(tibble())
    }
  }else{

    fit <- qlmer(formula = formula, data = fit_data, ...)

    tidy(fit$result) %>%
      mutate(metric = metric_name, model = list(fit$result), model_id = model_id, 
             n_ts = n_distinct(fit_data$rarefyID),
             formula = formula,
             warning = ifelse(length(fit$warnings) == 0, NA, fit$warnings))
  }
}

get_fixed_slopes <- function(metric, model, model_id, model_data){
  formula <- paste("pairwise", "~", model_id)
  fit_trends <- emtrends(model, specs = as.formula(formula), var="year_scaled",  mode = "satterth", lmerTest.limit = 20105, data = model_data)

  return(fit_trends)
}
````

```r
metric_names <- unique(model_data$metric)
raw_metrics <- c("SES_FRic", "SES_FDiv", "SES_FEve", 
                 "CWM_diet_inv", "CWM_forstrat_ground", "CWM_forstrat_watbelowsurf", "FEve")
log_metrics <- c("FRic", "S", "FDiv")
logone_metrics <- metric_names[!metric_names %in% c(raw_metrics, log_metrics)]

raw_value_mods <- tibble(metric_name = raw_metrics, 
                         base = "value ~ year_scaled + (year_scaled|study_id/rarefyID)",
                         taxa = "value ~ year_scaled * taxa + (year_scaled|study_id/rarefyID)",
                         realm = "value ~ year_scaled * realm + (year_scaled|study_id/rarefyID)",
                         climate = "value ~ year_scaled * climate + (year_scaled|study_id/rarefyID)") %>%
  pivot_longer(cols = c(base, taxa, realm, climate), names_to = "model_id", values_to = "formula") %>%
  #don't fit the CWM_litter_size_min_n and climate model, not enough data
  filter(metric_name != "CWM_litter_size_min_n" | model_id != "climate")

log_mods <- tibble(metric_name = log_metrics, 
                         base = "log(value) ~ year_scaled + (year_scaled|study_id/rarefyID)",
                         taxa = "log(value) ~ year_scaled * taxa + (year_scaled|study_id/rarefyID)",
                         realm = "log(value) ~ year_scaled * realm + (year_scaled|study_id/rarefyID)",
                         climate = "log(value) ~ year_scaled * climate + (year_scaled|study_id/rarefyID)") %>%
  pivot_longer(cols = c(base, taxa, realm, climate), names_to = "model_id", values_to = "formula")

logone_mods <- tibble(metric_name = logone_metrics, 
                         base = "logvalue ~ year_scaled + (year_scaled|study_id/rarefyID)",
                         taxa = "logvalue ~ year_scaled * taxa + (year_scaled|study_id/rarefyID)",
                         realm = "logvalue ~ year_scaled * realm + (year_scaled|study_id/rarefyID)",
                         climate = "logvalue ~ year_scaled * climate + (year_scaled|study_id/rarefyID)") %>%
  pivot_longer(cols = c(base, taxa, realm, climate), names_to = "model_id", values_to = "formula")

mods <- bind_rows(raw_value_mods, log_mods, logone_mods) %>% 
  select(metric_name, formula, model_id)

#Fit models with base settings
mod_results <- pmap_dfr(mods %>% filter(model_id == "base"), get_model, data = amph_data) %>%
  mutate(REML = TRUE, bobyqa = FALSE)

#manuscript table
mod_results %>% 
  filter(model_id == "base") %>% 
  select(-c(model_id, warning, statistic, df)) %>% 
  readr::write_csv(here::here("figures", "amph_model_table.csv"))