Model of Grampians Occupancy and Trait Data with linear relationships

knitr::opts_chunk$set(cache = TRUE, dev = c("png", "pdf"))
library(blme)
library(digest)
library(eucs)
library(here)
library(openblasctl)
library(xtable)

Data

mdg <- transform(
  modeldata_grampians,
  y   = occupancy,
  mlq = scale.default(log(mean_moisture_index_of_low_qtr)),
  twi = scale.default(log(topographic_wetness_index_3sec)),
  r1k = scale.default(log(relief_1000m_radius)),
  tn  = scale.default(log(total_nitrogen)),
  sla = scale.default(log(sla_mm2_per_mg)),
  sm  = scale.default(log(seed_mass_mg)),
  mh  = scale.default(log(max_height_m))
)

Model

call <- list(
  formula =
    y ~ mlq + twi + r1k + tn +
      mlq:sla + mlq:sm + mlq:mh +
      twi:sla + twi:sm + twi:mh +
      r1k:sla + r1k:sm + r1k:mh +
      tn:sla + tn:sm + tn:mh +
      (mlq + twi + r1k + tn | taxon),
  family = quote(binomial),
  data = mdg,
  fixef.prior = quote(normal(sd = 1)),
  cov.prior = quote(
    invwishart(
      df = 5 + 4, scale = diag(sqrt((5 + 4) / 2), 5), common.scale = FALSE
    )
  ),
  control = quote(
    glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e6))
  )
)
file <- paste0(eucs_cache_path(), "/", digest(deparse(call)), ".rds")

call$data <- quote(mdg)

if (!file_test("-f", file)) {
  openblas_set_num_threads(30)
  saveRDS(do.call(bglmer, call), file)
  openblas_set_num_threads(1)
}

model_grampians <- readRDS(file)

summary(model_grampians)

Cov prior  : taxon ~ invwishart(df = 9, scale = c(2.1213, 0, 0, ...), posterior.scale = cov, common.scale = TRUE)
Fixef prior: normal(sd = c(1, ...), corr = c(0 ...), common.scale = FALSE)normal(sd = c(1, ...), corr = c(0 ...), common.scale = FALSE)
Prior dev  : 25.4848

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [bglmerMod]
 Family: binomial  ( logit )
Formula: y ~ mlq + twi + r1k + tn + mlq:sla + mlq:sm + mlq:mh + twi:sla +  
    twi:sm + twi:mh + r1k:sla + r1k:sm + r1k:mh + tn:sla + tn:sm +  
    tn:mh + (mlq + twi + r1k + tn | taxon)
   Data: mdg
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))

     AIC      BIC   logLik deviance df.resid 
  4454.0   4681.8  -2195.0   4390.0     9108 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0649 -0.3144 -0.1498 -0.0643 15.0957 

Random effects:
 Groups Name        Variance Std.Dev. Corr                   
 taxon  (Intercept) 1.3668   1.1691                          
        mlq         0.4941   0.7029   -0.21                  
        twi         0.1453   0.3812   -0.02  0.43            
        r1k         0.2271   0.4765   -0.12 -0.36 -0.22      
        tn          0.1329   0.3645    0.05 -0.43 -0.26  0.16
Number of obs: 9140, groups:  taxon, 20

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.24440    0.26783 -12.114  < 2e-16 ***
mlq         -0.25659    0.18456  -1.390  0.16444    
twi          0.11136    0.10897   1.022  0.30681    
r1k          0.39596    0.14860   2.665  0.00771 ** 
tn          -0.10482    0.12962  -0.809  0.41869    
mlq:sla     -0.01301    0.35069  -0.037  0.97042    
mlq:sm       0.20270    0.31765   0.638  0.52340    
mlq:mh       0.05000    0.22033   0.227  0.82048    
twi:sla      0.41885    0.20805   2.013  0.04409 *  
twi:sm      -0.10689    0.18786  -0.569  0.56939    
twi:mh       0.03335    0.13206   0.253  0.80062    
r1k:sla      0.47362    0.26187   1.809  0.07052 .  
r1k:sm       0.62574    0.23248   2.692  0.00711 ** 
r1k:mh      -0.03077    0.16454  -0.187  0.85164    
tn:sla      -0.24740    0.24130  -1.025  0.30523    
tn:sm        0.09289    0.20842   0.446  0.65581    
tn:mh       -0.10817    0.14776  -0.732  0.46414    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Table of random effects.

covariate_nms <- c(
  "MMI",
  "TWI",
  "R1k",
  "TN"
)

vc <- summary(model_grampians)$varcor$taxon
stdev <- as.character(round(sqrt(diag(vc)), 1))
corr <- matrix(as.character(round(cov2cor(vc), 1)), 5)
corr[!lower.tri(corr)] <- ""
vtab <- cbind(stdev, corr)
vtab <- rbind(
  c("$\\sigma$", "$\\rho$", "", "", "", ""),
  rbind(c("", "$B_0$", covariate_nms)),
  vtab
)
vtab <- cbind(
  rownames(vtab) <- c("", "", "$B_0$", covariate_nms), vtab[, -6]
)

vtab <- xtable(
  vtab, caption = "Random effects", label = "tab:varcovlin",
  align = "lll|llll", digits = 1
)

dir.create(
  here("notebooks/model_linear_grampians_files/figure-html"),
  showWarnings = FALSE, recursive = TRUE
)

print(
  vtab,
  file = here(
    "notebooks/model_linear_grampians_files/figure-html/vtab.tex"
  ),
  include.rownames = FALSE,
  table.placement = "htbp!",
  caption.placement = "top",
  hline.after = c(-1, 2, nrow(vtab)),
  include.colnames = FALSE,
  sanitize.text.function = identity
)

Table of fixed effects.

fe <- round(summary(model_grampians)$coefficients[, 1:2], 1)
rownames(fe) <- c(
  "$B_0$",
  covariate_nms,
  paste(rep(covariate_nms, each = 3), c("SLA", "SM", "MH"), sep = "-")
)
colnames(fe) <- c("$\\mu$", "$\\sigma$")

fe <- xtable(
  fe, caption = "Fixed effects", label = "tab:felin",
  align = "lll", digits = 1
)

print(
  fe,
  file = here(
    "notebooks/model_linear_grampians_files/figure-html/fetab.tex"
  ),
  include.rownames = TRUE,
  table.placement = "htbp!",
  caption.placement = "top",
  hline.after = c(-1, 1, 5, nrow(fe)),
  include.colnames = TRUE,
  sanitize.text.function = identity
)