I have several large GAM models to fit (lots of data, various distributions for residuals and random effect smooths). Currently, I use foreach with a cluster cpus based on the parallel package to run each of the models using the gam function on an individual node. This is an improvement, but I am wondering if it is possible to use parallelization within the bam function itself while at the same time running each model using foreach. Ideally, I would like the results of the model runs to be as similar as possible to those of the original gam function. Should I use the cluster argument to bam or the nthreads argument?