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

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) +
      (1 | locality),
  family = quote(binomial),
  data = mdg,
  fixef.prior = quote(normal(sd = 1)),
  cov.prior = quote(
    list(
      taxon ~ invwishart(
        df = 5 + 4, scale = diag(sqrt((5 + 4) / 2), 5), common.scale = FALSE
      ),
      locality ~ invgamma
    )
  ),
  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  : locality ~ gamma(shape = 0.5, scale = 100, posterior.scale = sd, common.scale = TRUE)
           : 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  : 275.7107

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) + (1 | locality)
   Data: mdg
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))

     AIC      BIC   logLik deviance df.resid 
  4578.5   4813.5  -2256.3   4512.5     9107 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6494 -0.3130 -0.1386 -0.0506 21.6972 

Random effects:
 Groups   Name        Variance Std.Dev. Corr                   
 locality (Intercept) 0.6633   0.8144                          
 taxon    (Intercept) 1.5539   1.2466                          
          mlq         0.7327   0.8560   -0.16                  
          twi         0.1633   0.4041    0.00  0.47            
          r1k         0.2875   0.5362   -0.07 -0.48 -0.30      
          tn          0.1541   0.3925    0.00 -0.49 -0.34  0.27
Number of obs: 9140, groups:  locality, 457; taxon, 20

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.51997    0.28872 -12.192  < 2e-16 ***
mlq         -0.37441    0.22784  -1.643  0.10031    
twi          0.14201    0.12387   1.146  0.25160    
r1k          0.46892    0.17221   2.723  0.00647 ** 
tn          -0.11293    0.15567  -0.725  0.46818    
mlq:sla     -0.08339    0.40390  -0.206  0.83644    
mlq:sm       0.23857    0.36938   0.646  0.51837    
mlq:mh       0.03660    0.26758   0.137  0.89121    
twi:sla      0.44036    0.21916   2.009  0.04451 *  
twi:sm      -0.15039    0.19943  -0.754  0.45080    
twi:mh       0.03121    0.14060   0.222  0.82433    
r1k:sla      0.53850    0.28882   1.864  0.06225 .  
r1k:sm       0.67463    0.26002   2.595  0.00947 ** 
r1k:mh      -0.04656    0.18383  -0.253  0.80005    
tn:sla      -0.27857    0.25553  -1.090  0.27565    
tn:sm        0.13768    0.22270   0.618  0.53641    
tn:mh       -0.11256    0.15897  -0.708  0.47893    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
IycgLS0tCiMnIHRpdGxlOiAiTW9kZWwgb2YgR3JhbXBpYW5zIE9jY3VwYW5jeSBhbmQgVHJhaXQgRGF0YSB3aXRoIGxpbmVhciByZWxhdGlvbnNoaXBzIGFuZCBsb2NhdGlvbiByYW5kb20gZWZmZWN0IgojJyBhdXRob3I6ICJXaWxsaWFtIEsuIE1vcnJpcyIKIycgZGF0ZTogImByIFN5cy5EYXRlKClgIgojJyBvdXRwdXQ6CiMnICAgcm1hcmtkb3duOjpodG1sX25vdGVib29rOgojJyAgICAgY29kZV9mb2xkaW5nOiBoaWRlCiMnIC0tLQoKIysgc2V0dXAsIG1lc3NhZ2U9RkFMU0UKa25pdHI6Om9wdHNfY2h1bmskc2V0KGNhY2hlID0gVFJVRSwgZGV2ID0gYygicG5nIiwgInBkZiIpKQpsaWJyYXJ5KGJsbWUpCmxpYnJhcnkoZGlnZXN0KQpsaWJyYXJ5KGV1Y3MpCmxpYnJhcnkob3BlbmJsYXNjdGwpCgojJyAjIyBEYXRhCiMrIGxvYWQtZGF0YQptZGcgPC0gdHJhbnNmb3JtKAogIG1vZGVsZGF0YV9ncmFtcGlhbnMsCiAgeSAgID0gb2NjdXBhbmN5LAogIG1scSA9IHNjYWxlLmRlZmF1bHQobG9nKG1lYW5fbW9pc3R1cmVfaW5kZXhfb2ZfbG93X3F0cikpLAogIHR3aSA9IHNjYWxlLmRlZmF1bHQobG9nKHRvcG9ncmFwaGljX3dldG5lc3NfaW5kZXhfM3NlYykpLAogIHIxayA9IHNjYWxlLmRlZmF1bHQobG9nKHJlbGllZl8xMDAwbV9yYWRpdXMpKSwKICB0biAgPSBzY2FsZS5kZWZhdWx0KGxvZyh0b3RhbF9uaXRyb2dlbikpLAogIHNsYSA9IHNjYWxlLmRlZmF1bHQobG9nKHNsYV9tbTJfcGVyX21nKSksCiAgc20gID0gc2NhbGUuZGVmYXVsdChsb2coc2VlZF9tYXNzX21nKSksCiAgbWggID0gc2NhbGUuZGVmYXVsdChsb2cobWF4X2hlaWdodF9tKSkKKQoKIycgIyMgTW9kZWwKIysgcnVuLW1vZGVsLCBtZXNzYWdlPUZBTFNFLCBjYWNoZT1GQUxTRQpjYWxsIDwtIGxpc3QoCiAgZm9ybXVsYSA9CiAgICB5IH4gbWxxICsgdHdpICsgcjFrICsgdG4gKwogICAgICBtbHE6c2xhICsgbWxxOnNtICsgbWxxOm1oICsKICAgICAgdHdpOnNsYSArIHR3aTpzbSArIHR3aTptaCArCiAgICAgIHIxazpzbGEgKyByMWs6c20gKyByMWs6bWggKwogICAgICB0bjpzbGEgKyB0bjpzbSArIHRuOm1oICsKICAgICAgKG1scSArIHR3aSArIHIxayArIHRuIHwgdGF4b24pICsKICAgICAgKDEgfCBsb2NhbGl0eSksCiAgZmFtaWx5ID0gcXVvdGUoYmlub21pYWwpLAogIGRhdGEgPSBtZGcsCiAgZml4ZWYucHJpb3IgPSBxdW90ZShub3JtYWwoc2QgPSAxKSksCiAgY292LnByaW9yID0gcXVvdGUoCiAgICBsaXN0KAogICAgICB0YXhvbiB+IGludndpc2hhcnQoCiAgICAgICAgZGYgPSA1ICsgNCwgc2NhbGUgPSBkaWFnKHNxcnQoKDUgKyA0KSAvIDIpLCA1KSwgY29tbW9uLnNjYWxlID0gRkFMU0UKICAgICAgKSwKICAgICAgbG9jYWxpdHkgfiBpbnZnYW1tYQogICAgKQogICksCiAgY29udHJvbCA9IHF1b3RlKAogICAgZ2xtZXJDb250cm9sKG9wdGltaXplciA9ICJib2J5cWEiLCBvcHRDdHJsID0gbGlzdChtYXhmdW4gPSAxZTYpKQogICkKKQpmaWxlIDwtIHBhc3RlMChldWNzX2NhY2hlX3BhdGgoKSwgIi8iLCBkaWdlc3QoZGVwYXJzZShjYWxsKSksICIucmRzIikKCmNhbGwkZGF0YSA8LSBxdW90ZShtZGcpCgppZiAoIWZpbGVfdGVzdCgiLWYiLCBmaWxlKSkgewogIG9wZW5ibGFzX3NldF9udW1fdGhyZWFkcygzMCkKICBzYXZlUkRTKGRvLmNhbGwoYmdsbWVyLCBjYWxsKSwgZmlsZSkKICBvcGVuYmxhc19zZXRfbnVtX3RocmVhZHMoMSkKfQoKbW9kZWxfZ3JhbXBpYW5zIDwtIHJlYWRSRFMoZmlsZSkKCnN1bW1hcnkobW9kZWxfZ3JhbXBpYW5zKQo=