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=