# Choosing error-variance priors in hierarchical models

I am trying to find a reasonable and largely uninformative set of priors for the error variance in a multi-level model. The model was developed by others, and I am unsure whether their choices were optimal.

The data consist of $J$ observations for each of $N$ individuals (a typical data structure would be $J=8$ and $N=500$). The authors model the outcome $Y$ for individual $i \in \{1 ... N\}$ and object $j \in \{1 ... J\}$ as: $$Y_{ij} \sim \mathcal{N}(\alpha_i + \beta_i\theta_j,\tau_{ij}^{-1}),$$ where $\alpha_i$, $\beta_i$ and $\theta_j$ are parameters of interest, and $\tau_{ij}$ is a precision term (inverse of the variance), constructed as $\tau_{ij} =\tau_{i}\tau_{j}$. (The use of precision is probably due to the authors' use of JAGS to run the model -- in my case, it may be more useful to define the standard deviation, as in Stan.)

The authors further place gamma priors on the precision terms: $\tau_{j} \sim \Gamma(0.1,0.1)$, and $\tau_{i} \sim \Gamma(\nu,\omega)$. For the latter, individual-specific precision terms, gamma hyperpriors are also placed on the shape and rate parameters: $\nu \sim \Gamma(0.1,0.1)$, and $\omega \sim \Gamma(0.1,0.1)$.

My question is whether these priors are a good choice, and what would be a good alternative? Here are a few considerations:

1. The use of inverse gamma priors on scale parameters has received criticism (e.g. Gelman et al. BDA3, p. 130-132), on the basis that they are not truly uninformative, and the posterior distribution may be sensitive to the choice of shape and scale parameters. Gelman now recommends uniform, half-normal, or half-cauchy distributions.
2. In practice, the current model does seem to give good parameter estimates (judging by correlations with external (highly valid) data, and simulations on fake data).
3. While all other aspects of the model converge quickly (i.e. giving low Rhats), the variance terms ($\tau_{i}$ and $\tau_{j}$) seem to take much longer to converge (if at all).

I have several related questions:

1. Do the priors outlined above make sense? Will the use of hyperpriors mitigate the issue of sensitity to choice of shape and scale parameters for the gamma distribution?
2. If I wanted to largely keep the structure of priors (the variance still being a product of $i$- and $j$-level terms), but use half-cauchy distributions, what would be a good way to do so?
3. How complicated should the prior structure be, i.e. is it necessary/desirable with hyperpriors, or might a simpler option be sufficient? (I realize that this may be context dependent.)

More generally, I am looking for a reasonable and largely uninformative set of priors, and would be grateful for any advice (I am obviously new to this).

An alternative that works well in Stan is described in a vignette for the rstanarm package. In short, a prior on a KxK covariance matrix is repeatedly decomposed as

• Correlation matrix
• Cholesky factor of the correlation matrix
• Vector of variances
• Simplex vector for the k-th proportion of variance for all k
• Trace of covariance matrix
• K
• Square of a scale parameter

It is possible to put what Stan calls the LKJ prior on the Cholesky factor of a correlation matrix such that all correlation matrices of order K are equally probably a priori, which is maximally uninformative.

It is also possible to put a Dirichlet prior on the simplex vector that is uniform among simplex vectors of order K by setting all of the shape parameters equal to 1, which is maximally uniformative.

The variances are the product of the simplex vector and the scalar trace of the covariance matrix, which in turn is the product of K and the square of an unknown scale parameter. The rstanarm package actually does a Gamma prior on this scale parameter, but the Jeffreys prior would be considered more uniformative.

This combination typically results in good convergence, good effective sample size, and the posterior distribution is driven almost entirely by the data.

• Thank you, this seems like an excellent answer. It may be a little over my head, and I am struggling a bit to see how I would implement this. The examples I see in the rstanarm vignette and the Stan manual consider coefficients that are assumed to be unbounded and MVN distributed, while I want to model SDs with a lower bound at 0. (Note that I am interested in modeling heteroskedasticity and not autocorrelation.) Would something like a half MVN distribution (with a JxJ covariance matrix) be an option to model sigma_ij? (Apologies if I am a bit slow.) – notfound Feb 13 '17 at 13:23
• You could use these ideas for other things besides coefficients. It is probably easiest to use a MVN distribution to model the logarithm of each sigma_ij. – Ben Goodrich Feb 26 '17 at 4:05