Random Poisson variable is a negative but must be nonnegative!

You are expecting sampling statements like nJ[i,j] ~ poisson(a[i]*J[j]/totalArea); to populate nJ[i,j] with a value sampled from a Poisson distribution. But that’s not what these statements achieve in Stan. Instead, in Stan, these statements exist to define something proportional to the joint log-density (or log-mass) of the data and the parameters. I’m not too clear on why nJ[i,j] is getting initialized as -2147483648, but what you really need to do for latent parameters is to declare them as parameters, so that they will get initialized to some appropriate value subject to the appropriate constraints, and proceed from there.

However, this brings up a separate problem, because you cannot have discrete parameters in Stan. So to fit this model, what you really need is to marginalize nA and nJ out of the likelihood. Fortunately, these poisson-binomials marginalize straightforwardly to poissons, so you can just write, for example C[i,j,1] ~ poisson(c[i] * a[i] * J[j] / totalArea)

Note that this marginalization works only because each sampled poisson variate is only used once in the binomial likelihood. In cases where a Poisson sample serves as the number of counts in multiple binomial samples (e.g. in N-mixture models), you need an alternative scheme marginalize over the variate, e.g. by brute-force summing through the terms up to some bound of vanishingly low upper-tail probability in the Poisson.