Simplifying some code chunks in the Priors vignette

[strengejacke]

Some of the code chunks in the priors vignette could be simplified, in particular those related to create the plots.

WARNING:

1. the suggestions are heavily biased towards using easystats-framework / parameters-package
2. this means that the proposed changes would require adding see and parameters to the suggested packages in the description file.

library(glmmTMB)
library(lme4)
#> Loading required package: Matrix
library(blme)
library(ggplot2); theme_set(theme_bw())
library(parameters)
OkIt <- unname(palette.colors(n = 8, palette = "Okabe-Ito"))[-1]

cdat <- readRDS(system.file("vignette_data", "culcita.rds", package = "glmmTMB"))
cdatx <- cdat[-20, ]

form <- predation~ttt + (1 | block)
cmod_glmer <- glmer(form, data = cdatx, family = binomial)
cmod_glmmTMB <- glmmTMB(form, data = cdatx, family = binomial)
cmod_bglmer <- bglmer(form,
  data = cdatx, family = binomial,
  fixef.prior = normal(cov = diag(9, 4))
)

cprior <- data.frame(
  prior = rep("normal(0,3)", 2),
  class = rep("beta", 2),
  coef = c("(Intercept)", "")
)
print(cprior)
#>         prior class        coef
#> 1 normal(0,3)  beta (Intercept)
#> 2 normal(0,3)  beta

cmod_glmmTMB_p <- update(cmod_glmmTMB, priors = cprior)

all.equal(coef(summary(cmod_bglmer)),
  coef(summary(cmod_glmmTMB_p))$cond,
  tolerance = 5e-2
)
#> [1] TRUE

compare_parameters(
  cmod_glmer,
  cmod_glmmTMB,
  cmod_bglmer,
  cmod_glmmTMB_p
) |> plot(show_intercept = TRUE)

# errors
# confint(cmod_bglmer, method = "profile")


gdat <- readRDS(system.file("vignette_data", "gophertortoise.rds", package = "glmmTMB"))
form <- shells~prev + offset(log(Area)) + factor(year) + (1 | Site)
gmod_glmer <- glmer(form, family = poisson, data = gdat)
#> boundary (singular) fit: see help('isSingular')

gmod_bglmer <- bglmer(form, family=poisson, data=gdat)
## cov.prior = gamma(shape = 2.5, rate = 0, common.scale = TRUE, posterior.scale = "sd"))
gmod_glmmTMB <- glmmTMB(form, family = poisson, data = gdat) ## 1e-5
## bglmer default corresponds to gamma(Inf, 2.5)
gprior <- data.frame(
  prior = "gamma(1e8, 2.5)",
  class = "ranef"
)
gmod_glmmTMB_p <- update(gmod_glmmTMB, priors = gprior)
vc1 <- c(VarCorr(gmod_glmmTMB_p)$cond$Site)
vc2 <- c(VarCorr(gmod_bglmer)$Site)
all.equal(vc1, vc2, tolerance = 5e-4)
#> [1] TRUE

# profiled CI error for gmod_bglmer
compare_parameters(
  gmod_glmer, gmod_glmmTMB, gmod_bglmer, gmod_glmmTMB_p,
  ci_method = "profile",
  effects = "random"
) |> plot(show_intercept = TRUE)

# works w/o error
compare_parameters(
  gmod_glmer, gmod_glmmTMB, gmod_bglmer, gmod_glmmTMB_p,
  effects = "random"
) |> plot(show_intercept = TRUE)

^{Created on 2024-04-05 with reprex v2.1.0}

Notes:

• I would add a sentence or two at one point or another to make it read less technically.
• Unlike you wrote in the vignette, profiled CI seems to work for glmer, but not for bglmer for me?
• I would update the data frame with the prior definition, because "theta" is not mentioned in the docs (see my example above).

WDYT? I could open a PR if you're interested in this change, but I also have no problem if you don't think this is necessary or don't want to add further "dependencies" to your package.