The big problem, as @betanalpha is pointing out, is that the log isn’t an adequate link function for probabilities. So what happens is that you chop out bits of the probability space that don’t lead to values in (0, 1). This might not be a problem in practice—sometimes you can just do simple linear regressions for probabilities if everything’s constrained enough. If it does crop up, what happens is that Stan has to do a rejection when it steps into an illegal region. Also, you have to be careful with normalization if the volume of the illegal region depends on parameters.
The compound declare define makes things cleaner:
int Kc = K - 1;
...
vector[N_1] r_1_1 = sd_1[1] * z_1[1];
You can turn th definition of `mu` into a one-liner by adding elementwise products and multi-indexing.
vector[N] mu = Xc * b + temp_Intercept + r_1_1[J_1] .* Z_1_1;
Ther emight be a way to remove some redundant multiplies in that `r_1_1` and `Z_1_1` term depending on how `J_1` is structured. The point is to eliminate redundant/repetitive computations whenever possible.
I'll leave the two expressions in teh generated quantities up to you.
I'm not sure why you only give one entry of `z_1` a prior here:
z_1[1] ~ normal(0, 1);
The easiest way to deal with "prior only" conditions is to just have size zero data.