some comments on your comments :):
I completely agree, in my opinion the implementation that you are talking about is perfectly right.
I also agree. Setting the “precision” to 1 does not mean to set the standard deviation to 1 in ICAR random effects, in part because each value of the random effect has a different standard deviation. But this also happens on the “precision” parameter tau as this parameter is multiplied in the diagonal of the precision matrix by the number of neighbours so we do not have a single common precision for all the observations. This is why I have quoted the term “precision” at the beginning of this paragraph. Anyway, when I state setting the precision or standard deviation I mean fixing the scale of the ICAR distribution what means in some sense to fix either the “precision” or the “standard deviation” (in a broad sense) to a certain value.
As mentioned below the original implementation did not include the sigma_phi^2 in the target term so the model is not doing exactly that. Anyway we thought about the modelling alternative that you are proposing. Just one comment about that, in the expression for log(p(phi)) that you put above there is one term left corresponding to the (determinant of the precision matrix)^(-1/2). The log-determinant of sigma^2(D-W) can be put as log(sigma)2(N-1) times the log-determinant of D-W (* I put two technical asides on this at the next paragraph). The second of these terms is just a constant but the first of them do matters for the log-likelihood. Thus implementing this proposal in terms of a “standard deviation” sigma_phi, in our opinion, should be done as:
target += -0.5 * dot_self((sigma_phi^2)*(phi[node1] - phi[node2])) - (N-1) *log(sigma_phi);
About the technical asides: the N-1 in (N-1) *log(sigma_phi) is valid if the region of study has a single connected component. If it was composed of K isolated connected components the N-1 should be changed to N-K. Besides, D-W is singular so its log-determinant is (minus) infinite. Nevertheless, the sum-to-zero restriction leaves out the null dimension of that matrix and the determinant of D-W could be changed by a generalized determinant, which turns back a finite constant value that we will not have to take care of. That is why that term is nos included in the last expression.
Anyway, at the end, we did not implement this proposal since this would be a centered parameterization of the ICAR random effects in the BYM model, what seems to be an inadvisable option according to the Stan manual. So we opted by other different modeling options with explicit non-centered parameterizations.
Ok, I agree that this could be improved. We are on that. Regretfully INLA is not an option for many of the models that I am handling, this is why we are trying Stan and considering some other MCMC options. Anyway, for now, what amazes us is the competitiveness of WinBUGS which I expected to be beated by far for any alternative implementation of the BYM model, for example, with more modern MCMC tools. For now, with our still preliminary results, WinBUGS unexpectedly keeps being a competitive option. Anyway, as you mention, I am also sure that more efficient Stan implementations are possible and Stan will be able to show all its power also in disease mapping problems. We are on that ;).
Thanks for your answers and your time :).