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

knitr::opts_chunk$set(echo = TRUE)

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

Get data and put it in long form

pins::board_register_github(repo = "karinorman/biodivTS_data", branch = "master")

model_data <- pins::pin_get("model-data", board = "github")

Data exploration

Visualize distributions of response

# model_data %>%
#   filter(metric %in% c("FRic", "FEve", "FDiv", "FDis", "S", "Jaccard_base")) %>%
#   ggplot() +
#   geom_histogram(aes(x = logvalue),
#                  bins = 60) +
#   facet_wrap(~metric, scales = "free") +
#   labs(title = "Distribution of Log values")
# 
# model_data %>%
#   filter(metric %in% c("FRic", "FEve", "FDiv", "FDis", "S", "Jaccard_base")) %>%
#   ggplot() +
#   geom_histogram(aes(x = value),
#                  bins = 60) +
#   facet_wrap(~metric, scales = "free") +
#   labs(title = "Distribution of raw values")

Fit a single model to play around with

test_data <- model_data %>% 
  filter(metric == "FRic") %>%
  mutate(year_scaled = scale(year),
           duration_scaled = scale(duration, scale = FALSE),
           startyear_scaled = scale(startyear, scale = FALSE),
         taxa = as.factor(taxa))

Different options for getting between group samples, fixed effects (test differences between intercept):

library(emmeans)

#fixed effects contrasts
fit_fixed <- lmer(formula = "value ~ year_scaled*taxa + (year_scaled|study_id/rarefyID)", data = test_data)

#intercept contrasts
m_emm <-  emmeans(fit_fixed,  mode = "satterth", lmerTest.limit = 20105, pairwise ~ taxa | year_scaled)
fit_cont <- as_tibble(m_emm$contrasts)

#slope contrasts and within group slope estimates
fit_trends <- emtrends(fit_fixed, pairwise ~ taxa, var="year_scaled",  mode = "satterth", lmerTest.limit = 20105) %>% as_tibble()

interactions::interact_plot(fit_fixed, pred = year_scaled, modx = taxa, plot.points = TRUE)

As a random effect:

fit_rand <- lmer(formula = "value ~ year_scaled + (year_scaled|study_id/rarefyID) + (year_scaled|taxa)", data = test_data)

#lognormal
ln_fit <- lmer(logvalue ~ year_scaled + (year_scaled|study_id/rarefyID), data = test_data)

gamma_fit <- glmer(value ~ year_center + (year_center|study_id/rarefyID), data = test_data, family = Gamma)

poiss_fit <- glmer(value ~ year_center + (year_center|study_id/rarefyID), data = test_data, family = poisson(link = "log"))

nb_fit <- glmer.nb(value ~ year_center + (year_center|study_id/rarefyID), data = test_data)

Get p-values for blups

cV <- ranef(ln_fit, effects = "ran_vals", drop = TRUE)
ranvar <- attr(cV[[1]], "postVar")

ng <- dim(ranvar)[3]
np <- dim(ranvar)[2]
mm <- matrix(ranvar[cbind(rep(seq(np),ng),
             rep(seq(np),ng),
             rep(ng,each=np))],
       byrow=TRUE,
       nrow=ng)

sumVar <- vcov(ln_fit)[1,1]+mm[,1]

year_scaled_var <- tidy(ln_fit) %>%
  filter(term == "year_scaled") %>% 
  mutate(var = std.error*std.error) %>%
  pull(var) 

study_ests <- broom.mixed::tidy(ln_fit, effects="ran_vals") %>%
  filter(group == "study_id", term == "year_scaled") %>% 
  rename(cond.std.error = std.error) %>%
  mutate(cond.var = cond.std.error*cond.std.error, 
         var = cond.var + year_scaled_var,
         std.error = sqrt(var),
         upr.ci = estimate + (2*std.error),
         lwr.ci = estimate - (2*std.error),
         sig = case_when(
           lwr.ci < 0 & upr.ci > 0 ~ FALSE,
           TRUE ~ TRUE
         ))

Look for best optimizer using allFit(), doesn't look like any of them work

#find best optimizer
#lognormal fit
ln_optims <- allFit(ln_fit, maxfun = 1e5)

is.OK <- sapply(ln_optims, is, "merMod")
ln_optims.OK <- ln_optims[is.OK]
lapply(ln_optims.OK,function(x) x@optinfo$conv$lme4$messages)

Let's try all the options for both optimix and nloptwrap, first for the generalized option w/Gamma distribution

optimx_options <- c("L-BFGS-B", "nlminb", "nlm", "bobyqa", "nmkb", "hjkb")

for(i in 1:length(optimx_options)){
  model_flex <- lmer(logvalue ~ year_center + (year_center|study_id/rarefyID), data = test_data,
                     control = lmerControl(optimizer = "optimx",
                                           optCtrl = list(method = optimx_options[i],
                                                                   maxit = 1e9)))
  if(is.null(model_flex@optinfo$conv$lme4$messages)){
    print(paste0("One of the optimx options, ", optimx_options[i],", worked!"))
    print(summary(model_flex))
    break
  }
}

algoptions <- c("NLOPT_LN_PRAXIS", "NLOPT_GN_CRS2_LM",
"NLOPT_LN_COBYLA", "NLOPT_LN_NEWUOA",
"NLOPT_LN_NEWUOA_BOUND", "NLOPT_LN_NELDERMEAD",
"NLOPT_LN_SBPLX", "NLOPT_LN_BOBYQA")

for(i in 1:length(algoptions)){
  model_flex <- lmer(logvalue ~ year_center + (year_center|study_id/rarefyID), data = test_data,
                     control = lmerControl(optimizer = "nloptwrap",
                                           optCtrl = list(algorithm = algoptions[i],
                                                          maxeval = 1e7,
                                                          xtol_abs = 1e-9,
                                                          ftol_abs = 1e-9)))
  if(is.null(model_flex@optinfo$conv$lme4$messages)){
    print(paste0("One of the nloptwrap options, ", algoptions[i],", worked!"))
    print(summary(model_flex))
    break
  }
}

optimx_options <- c("L-BFGS-B", "nlminb", "nlm", "bobyqa", "nmkb", "hjkb")

for(i in 1:length(optimx_options)){
  model_flex <- glmer(value ~ (year|study_id/rarefyID), data = test_data, family = Gamma,
                     control = glmerControl(optimizer = "optimx",
                                           optCtrl = list(method = optimx_options[i],
                                                                   maxit = 1e9)))
    if(is.null(model_flex@optinfo$conv$lme4$messages)){
    print(paste0("One of the optimx options, ", optimx_options[i],", worked!"))
    print(summary(model_flex))
    break
  }
}

#nloptwrap
algoptions <- c("NLOPT_LN_PRAXIS", "NLOPT_GN_CRS2_LM",
"NLOPT_LN_COBYLA", "NLOPT_LN_NEWUOA",
"NLOPT_LN_NEWUOA_BOUND", "NLOPT_LN_NELDERMEAD",
"NLOPT_LN_SBPLX", "NLOPT_LN_BOBYQA")

for(i in 1:length(algoptions)){
  model_flex <- glmer(value ~ (year|study_id/rarefyID), data = test_data, family = Gamma,
                     control = glmerControl(optimizer = "nloptwrap",
                                           optCtrl = list(algorithm = algoptions[i],
                                                          maxeval = 1e7,
                                                          xtol_abs = 1e-9,
                                                          ftol_abs = 1e-9)))
  if(is.null(model_flex@optinfo$conv$lme4$messages)){
    print(paste0("One of the nloptwrap options, ", algoptions[i],", worked!"))
    print(summary(model_flex))
    break
  }
}

test_data <- model_data %>% 
  filter(metric == "SES_FEve") %>%
  mutate(study_id = as.factor(study_id))

all_model <- lmer(value ~ year_scaled + taxa + realm + (year_scaled|study_id/rarefyID), 
                                 data = test_data)