I discovered tonight that I’ve spent years misunderstanding the example in the
SUG 1.13 (“Multivariate priors for hierarchical models”), and now have a couple questions:

The absence of a prior on
sigma
(observationlevel noise scale) is surely an error, yes? I can’t imagine that the implieduniform(0,+Inf)
could be expected to perform appropriately. 
I feel like there are serious identification issues in the model at present due to the combination of three structures:
a) The dataY
consist of single observations from each ofN
individuals, modelled asY~normal(mu,sigma);
b)mu
is itself a vector of lengthN
, each entry of which encodes the dotproduct between the corresponding row in the N \times K contrast matrixx
and the column in the K \times J coefficient matrixbeta
corresponding to the group of that individual.
c) The K \times J coefficient matrixbeta
is modelled as a multivariate normal with the square of the scale vectortau
forming the diagonal.
Thus, I don’t think that sigma
and tau
are actually identifiable with this structure. These quantities would be properly identified if we had multiple observations from each of the individuals (or even just one of the individuals), but that’s explicitly not how the data declarations are set up.
I actually think I long ago worked out (with help from other Stan gurus as I recall) code for a version of the model that doesn’t have this identification issue (old video walkthrough here, but the critical bit is the use of the get_coefs()
function in the transformed data section to transform the data to the coefficients that are then modelled as multivariate normal in the model).
Or am I missing something (again)?
(Tagging @andrewgelman but happy to receive input/correction from anyone)