# How can I run a Zero-Inflated Poisson/Negative Binomial Mixed Model with Gaussian Process

After having visited stats stack exchange countless times, I'm finally asking a question!

For my research, I am try to run a model of the form: $$Y = f(X,B)+ g(X) + \epsilon$$

Where $$f(X,B)$$ is a zero-inflated poisson/negative binomial density, $$g(X)$$ is a Gaussian process representing spatial correlation, and $$epsilon$$ is your usual error. Sorry if my notation is a little loose. Additionally, I want my zero-inflated portion to be able to have random effects. In summary, I need a program/package that can run:

1. Zero-inflated GLMs
2. Covariance structure (gaussian process)
3. Random Effects

However, I can only find programs that are capable of two out of these three items. R is highly preferred, but I am also open to SAS as a last resort. I already looked at looked at PROC GLMMIX in SAS, but that does not have all three of these items. I'm already using packages like pscl (which does not have all three items in the list above). Stan is too slow to handle the size of the dataset. I have around 100,000 observations, and it would appear that the processing time is on the order of days or even weeks (100 observations took 18 minutes, and this was with assigning each of the two chains I was running to its own core (i5). And this was before even using an random or spatial effects in the ZIP model! In my literature review, it would appear that R-INLA could potentially solve this issue if I stick with bayesian, but it would appear people still need computing clusters to make this work.

Can anyone suggest ways to resolve this ?

Background: I am modeling the abundance of a species of bird in Wisconsin. As you can imagine, most observations even for commons species have "zero" as the count. Hence the need for zero inflation. There is also observer data with repeated measures, hence the mixed model. And there is spatial correlation to deal with, hence the Gaussian process.