Model of Grampians Occupancy and Trait Data with quadratic 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(poly(log(mean_moisture_index_of_low_qtr), 2))[, 1],
  mlq2 = scale.default(poly(log(mean_moisture_index_of_low_qtr), 2))[, 2],
  twi  = scale.default(poly(log(topographic_wetness_index_3sec), 2))[, 1],
  twi2 = scale.default(poly(log(topographic_wetness_index_3sec), 2))[, 2],
  r1k  = scale.default(poly(log(relief_1000m_radius), 2))[, 1],
  r1k2 = scale.default(poly(log(relief_1000m_radius), 2))[, 2],
  tn   = scale.default(poly(log(total_nitrogen), 2))[, 1],
  tn2  = scale.default(poly(log(total_nitrogen), 2))[, 2],
  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 + mlq2 + twi + twi2 + r1k + r1k2 + tn + tn2 +
      mlq:sla + mlq:sm + mlq:mh + mlq2:sla + mlq2:sm + mlq2:mh +
      twi:sla + twi:sm + twi:mh + twi2:sla + twi2:sm + twi2:mh +
      r1k:sla + r1k:sm + r1k:mh + r1k2:sla + r1k2:sm + r1k2:mh +
      tn:sla + tn:sm + tn:mh + tn2:sla + tn2:sm + tn2:mh +
      (mlq + mlq2 + twi + twi2 + r1k + r1k2 + tn + tn2 | taxon),
  family = quote(binomial),
  data = mdg,
  fixef.prior = quote(normal(sd = 1)),
  cov.prior = quote(
    invwishart(
      df = 9 + 4, scale = diag(sqrt((9 + 4) / 2), 9), 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 = 13, scale = c(2.5495, 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  : -36.3444

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

     AIC      BIC   logLik deviance df.resid 
  4138.8   4694.2  -1991.4   3982.8     9062 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3002 -0.2788 -0.1186 -0.0389 15.6177 

Random effects:
 Groups Name        Variance Std.Dev. Corr                                     
 taxon  (Intercept) 1.16937  1.0814                                            
        mlq         0.21972  0.4687   -0.12                                    
        mlq2        0.56698  0.7530   -0.08 -0.47                              
        twi         0.11733  0.3425   -0.09  0.20  0.00                        
        twi2        0.08609  0.2934    0.13  0.03 -0.17 -0.07                  
        r1k         0.34319  0.5858   -0.29 -0.41  0.57 -0.18 -0.07            
        r1k2        0.14525  0.3811    0.39  0.17 -0.20  0.08  0.11 -0.18      
        tn          0.15571  0.3946   -0.08 -0.43  0.45 -0.17 -0.10  0.38 -0.27
        tn2         0.11103  0.3332   -0.08  0.04 -0.03  0.19 -0.04 -0.22 -0.13
      
      
      
      
      
      
      
      
      
 -0.09
Number of obs: 9140, groups:  taxon, 20

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.952517   0.270563 -14.609  < 2e-16 ***
mlq         -0.754599   0.238553  -3.163 0.001560 ** 
mlq2        -0.804486   0.223215  -3.604 0.000313 ***
twi          0.248615   0.124941   1.990 0.046607 *  
twi2        -0.230252   0.097996  -2.350 0.018793 *  
r1k          0.331276   0.244286   1.356 0.175068    
r1k2        -0.513143   0.156213  -3.285 0.001020 ** 
tn           0.038942   0.173418   0.225 0.822326    
tn2         -0.232909   0.116129  -2.006 0.044898 *  
mlq:sla     -0.648563   0.386835  -1.677 0.093623 .  
mlq:sm       0.234553   0.334256   0.702 0.482856    
mlq:mh       0.199616   0.249239   0.801 0.423188    
mlq2:sla    -0.624993   0.388772  -1.608 0.107921    
mlq2:sm     -0.314921   0.361699  -0.871 0.383934    
mlq2:mh     -0.169287   0.263254  -0.643 0.520187    
twi:sla      0.537087   0.218241   2.461 0.013856 *  
twi:sm      -0.062108   0.198786  -0.312 0.754708    
twi:mh      -0.015197   0.141515  -0.107 0.914483    
twi2:sla    -0.288002   0.178871  -1.610 0.107375    
twi2:sm     -0.120662   0.162668  -0.742 0.458226    
twi2:mh     -0.043196   0.113130  -0.382 0.702592    
r1k:sla      0.364892   0.322784   1.130 0.258287    
r1k:sm       0.877609   0.302484   2.901 0.003716 ** 
r1k:mh      -0.268953   0.212709  -1.264 0.206080    
r1k2:sla    -0.040378   0.225307  -0.179 0.857771    
r1k2:sm     -0.090706   0.205815  -0.441 0.659418    
r1k2:mh      0.099557   0.142949   0.696 0.486147    
tn:sla      -0.142778   0.278619  -0.512 0.608338    
tn:sm        0.090155   0.254761   0.354 0.723428    
tn:mh       -0.244277   0.182251  -1.340 0.180137    
tn2:sla     -0.000054   0.209273   0.000 0.999794    
tn2:sm       0.114641   0.190812   0.601 0.547969    
tn2:mh      -0.061827   0.129849  -0.476 0.633969    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Table of random effects.

covariate_nms <- c(
  mlq  = "MMI",
  mlq2 = "MMI$^2$",
  twi  = "TWI",
  twi2 = "TWI$^2$",
  r1k  = "R1k",
  r1k2 = "R1k$^2$",
  tn   = "TN",
  tn2  = "TN$^2$"
)

vc <- summary(model_grampians)$varcor$taxon
stdev <- as.character(round(sqrt(diag(vc)), 1))
corr <- matrix(as.character(round(cov2cor(vc), 1)), 9)
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[, -10]
)

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

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

print(
  vtab,
  file = here(
    "notebooks/model_quadratic_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:fequa",
  align = "lll", digits = 1
)

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