Gotchya. My understanding (which is probably wrong) is that arraydist
and filldist
sped things up when dealing with parameters, but not for the observations. Hopefully someone will correct me if I am wrong.
You could write it using standard broadcasting syntax. That worked fine for me when I rewrote your model 2.
θ = staffweights[staffid] .+ patientweights[patientid] .+ usewch .* wcwt
totweight .~ Gamma.(15, θ ./ 14)
Another option might be a custom distribution. An example using MvBinomial is here.
My only other (probably unhelpful) suggestions at this point to keep things positive are the obvious options (e.g., log(y) or using a something like a Gamma distribution, as in your example), which I assume you have tried or cannot use.