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.

I think that in this instance (if you are trying to fit a hurdle model count model) that you could use a single model in {mgcv} with the family = ziplss() option. Yes, from the name this family looks like it is a zero-inflated Poisson but it is not that zero-inflated Poisson. I really is a hurdle model - the Poisson part can't return a 0.
In that case, you'll need to specify a list of 2 formulas instead of one:
m <- gam(list(abund ~ s(VarR, k = 10),  # count part
                    ~ s(VarR, k = 10)), # binomial hurdle
         data = landcover, 
         method = "REML",               # <-- only REML!
         family = ziplss())

Note that with these general families that are way outside the common-or-garden exponential family of distributions that only REML smoothness selection is possible and hence you shouldn't be comparing models with different fixed effect parts (i.e. from your other question about AICc selection for these models.)
As I don't subscribe to that [your AICc approach in the other Q] kind of modelling approach, fit the model as per your prior expectations or knowledge from the literature and let it be estimated, or, you could add select = TRUE to add more penalties to the smooths such that entire smooth terms can be shrunk out of the model.
Also note the the caution in ?ziplss about excessive 0s in parts of the covariate space and what Simon Wood advises users to do if that is the case for their data.
If you are doing something more complicated like trying to include detection probabilities in your model for abundance, then consider the paper by Mark Bravington, Dave Miller, and Sharon Hedley (2021), which shows how to do this kind of thing with {mgcv}'s GAMs with the case study of distance sampling derived data.
References
Bravington, M.V., Miller, D.L., Hedley, S.L., 2021. Variance Propagation for Density Surface Models. J. Agric. Biol. Environ. Stat. https://doi.org/10.1007/s13253-021-00438-2
