Autoscale parameter tau in `decov()` prior?

Summary:

It might be necessary/better to autoscale the parameter "tau" in the decov() prior.

Description:

If I understand the documentation here, here, and here correctly, parameter tau (the key parameter for the trace of the covariance matrix of group-level terms) is not autoscaled in the decov() prior, in contrast to the corresponding parameter in the lkj() prior used by default in stan_mvmer() and stan_jm(). However, I think it would make sense to perform this autoscaling. Consider a linear regression with group-level intercepts. The default exponential prior for tau with rate parameter 1 (not autoscaled) would imply that a value of 3 for the standard deviation of the group-level intercepts would already be extreme (taking qexp(0.95)). However, the standard deviation of y (and correspondingly also the true standard deviation of the group-level intercepts) may be much larger than 3 (in fact, it may be arbitrarily large; for example, y may have a standard deviation of 50).

The brms package also performs that autoscaling, see e.g.

data("kidiq", package = "rstanarm")
kidiq_gr <- within(kidiq, {
  mom_age_gr <- cut(mom_age,
                    breaks = unique(quantile(mom_age, probs = seq(0, 1, 0.1))),
                    include.lowest = TRUE)
})
brms::make_stancode(formula = kid_score ~ mom_hs + mom_iq + (1 | mom_age_gr),
                    data = kidiq_gr)

which gives in the model block:

[...]
  // priors including constants
  target += student_t_lpdf(Intercept | 3, 90, 19.3);
  target += student_t_lpdf(sigma | 3, 0, 19.3)
    - 1 * student_t_lccdf(0 | 3, 0, 19.3);
  target += student_t_lpdf(sd_1 | 3, 0, 19.3)
    - 1 * student_t_lccdf(0 | 3, 0, 19.3);
[...]

where sd_1 is tau mentioned above and 19.3 is mad(kidiq$kid_score).

The documentation here states:

If all the variables were multiplied by a number, the trace of their covariance matrix would increase by that number squared. Thus, it is reasonable to use a scale-invariant prior distribution for the positive scale parameter, [...]

but I'm not sure if this is correct. Can we really arbitrarily multiply all the variables by a number? In case of only group-level intercepts for example, these variables are the "auxiliary" group-level intercepts (i.e., the group-level intercepts before multiplying them with their standard deviation) and are therefore assumed to follow a unit normal distribution. So they are restricted in their scale and it's the multiplication with tau which enables them to take on arbitrarily large values. I don't think there's an identifiability issue here.

I'm not too familiar with rstanarm's internals, but I think this autoscaling issue could be resolved by adding an autoscale argument to the decov() function. See also the reproducible steps below.

Reproducible Steps:

options(mc.cores = parallel::detectCores(logical = FALSE))
data("kidiq", package = "rstanarm")
kidiq_gr <- within(kidiq, {
  mom_age_gr <- cut(mom_age,
                    breaks = unique(quantile(mom_age, probs = seq(0, 1, 0.1))),
                    include.lowest = TRUE)
})
library(rstanarm)
myfit <- stan_glmer(kid_score ~ mom_hs + mom_iq + (1 | mom_age_gr),
                    data = kidiq_gr,
                    seed = 1466333131)
prior_summary(myfit)

gives as output:

[...]
Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
[...]

so the decov() prior indeed doesn't seem to be autoscaled which may also be seen from the fact that the following does not work:

myfit_autoscale <- update(
  myfit,
  prior_covariance = decov(regularization = 1, concentration = 1,
                           shape = 1, scale = 1,
                           autoscale = TRUE)
)

(as it throws the error

Error in decov(regularization = 1, concentration = 1, shape = 1, scale = 1,  : 
  unused argument (autoscale = TRUE)

).

RStanARM Version:

2.21.2 (from the Stan repo)

R Version:

4.1.0

Operating System:

Ubuntu 20.04.2 LTS

Start with the decov() prior documentation and its handling through stan_glmer(..., prior_covariance=...), then compare it with the autoscaling behavior described for lkj(). Reproduce the kidiq example and the autoscale argument error; done means the requested decov() autoscaling behavior is supported and reflected in the prior summary and documentation.