Thanks Martin and Bob!
Bob, I order those two parameters to identify the model because multiplying the parameters zeta and theta by -1 will result in the same relative distances between items (zeta and theta just represent spatial positions on a unidimensional scale) and I’ll get a bimodal posterior. I’m essentially arbitrarily choosing the direction of the scale, such that one item zeta[1] in the latent space that I’m trying to model is always larger than zeta[2]. I’ve tried identifying the model in alternative ways, such as fixing zeta[1] to +1 and zeta[2] to -1, and constraining the sign of one (or more) of the zeta. (http://www.stat.columbia.edu/~gelman/research/published/171.pdf, pp. 176-178) In any case, I’ve tried a variety of ways of doing this, and I still sometimes get a chain that explores the “wrong” part of the parameter space and doesn’t break out even with long warm-ups. With reasonable initializations, it all works as I suspect it should and the results make a lot of sense.
The .stan file is the centered version, but I have a non-centered version, which unfortunately still runs into the same problem.
Martin, I’m going to trying simulating the DGP to make sure I’ve coded this correctly. I usually do that when coding a model in Stan, but I skipped it this time, and perhaps I shouldn’t have. I’ll respond to this thread with the results.
Martin and Bob, thanks also for pointing out the ways to clean up the code. I’ll implement those changes to clean this up.