If I am not wrong, you are fitting a linear mixed model with two i.i.d. random effects (on intercept).
There is no need to use gamm
in this case. Using gam
with method = REML
will do.
gam(metric1~a+b+c+d+e+f+g+h+i+ s(X, bs = 're') + s(Y, bs = 're'), data = dataset, method = 'REML')
Note I have not extended other fixed effects into smooth functions, which you may do it yourself.
If you have a large dataset, it is suggested to use bam
, and note it is method = 'fREML
in this case.
The difference between gam
and bam
in REML estimation, is that the former uses "outer" iteration, while the latter uses "performance" iteration. But for Gaussian data there's no difference, though bam
itself is designed for large dataset, using iterative QR reduction and parallel computing on request.
Personally I think gamm
is outdated. It does REML estimation using lme
and MASS::glmmPQL
, which is much less efficient then the penalized least squares adopted by gam
and bam
.
mgcv::gam
– borgs