For looping lpmf’s is much slower than vectorized version. Both bernoulli_lpmf and neg_binomial_2_lpmf accept vector/array arguments, so there is no need to have for loop. See related example in thread Variable Selection in Parametric Survival Analysis Models
The same thread also discusses easy speedups by using appropriate stanc and C++ optimization flags. The same thread discusses also further speedups. As @paul.buerkner mentioend, swithcing to slightly simpler model would help, too.