How to multiply two models together in R? I am constructing a two-stage GAM to model capture rate across my study site while accounting for zero inflation.
VarRPA <- gam(p.a ~ s(VarR, k=10), data=landcover, method = "ML", 
              family = "binomial"(link = "logit"))
VarRAbund <- gam(abund ~ s(VarR, k=10), data=landcover, 
                 method = "ML", family = "poisson"(link = "log"))

For my overall prediction, I need to multiply the models by each other i.e. make them a single model rather than two subsequent ones. I know the hurdle function in the pscl package allows you to do this by effectively considering the hurdle model as a single model however I've constructed my models using the mgcv package and would, for several reasons, prefer to keep it that way. To get this overall prediction is it as simple as multiplying for each data point the outcome of stage 1 with the outcome of stage 2? I can do this manually in excel of course, but is there a way to get R to do it for me by reading it as a single model? This would also be incredibly helpful when it comes to cross-validating the model as so far I only know how to cross-validate each stage separately rather than how to cross-validate the multiplied model.

There are a variety of zero-inflated [ZI],
and zero-adjusted [ZA] (i.e. Hurdle) distributions available in the gamlss R package.
These include:
ZAP and ZIP, where P = Poisson
ZALG, where LG = logarithmic
ZAZIPF, where ZIPF = zipf
ZANBI and ZINBI, where NBI = negative binomial
ZAPIG and ZIPIG, where PIG = Poisson-inverse Gaussian
ZASICHEL and ZISICHEL, where SICHEL = Sichel,
ZABNB and ZIBNB, where BNB = beta negative binomial
With all the above there is no need to fit 2 separate models.
