In addition to mgcv and its zero-inflated Poisson families (ziP()
and ziplss()
), you might also look at the brms package by Paul-Christian Bürkner. It can fit distribution models (where you model more than just the mean, in your case the zero-inflation component of the model can be modelled as a function of covariates just like the count function).
You can include smooths in any of the linear predictors (for the mean/count, zero-inflation part, etc) via s()
and t2()
terms for simple 1-d or isotropic 2-d splines, or anisotropic tensor product splines respectively. It has support for zero-inflated binomial, Poisson, negative binomial, and beta distributions, plus zero-one-inflated beta distributions. It also has hurdle models for Poisson and negative binomial responses (where the count part of the model is a truncated distribution so as to not produce further zero counts).
brms fits these models using STAN, so they are fully Bayesian, but this will require you to learn a new set of interfaces to extract relevant information. That said, there are several packages offering support functions for just this task, and brms has helper functions written that utilise these secondary packages. You'll need to get STAN installed and you'll need a C++ compiler as brms compiles the model as defined using R into STAN code for evalutation.