GAMs: specifying knots positions for thin-plate regression splines PROBLEM STATEMENT
I have many replicates of the same dataset. I want to fit the same GAM to each replicate, and then average the model coefficients and covariance matrices (Vp) from all GAMs. The best smooth for my variables is a thin-plate regression spline. However, when running the GAMs with bs = "tp", there is not information available about the position of the knots for the splines in the model. When specifying the knots, the AIC is much higher than that of the GAM where the algorithm positions the knots (R-sq also drops substantially when knots positions are given).
QUESTION
Is it possible to have a GAM in which I specify the knots positions, but without having a big loss of model fit? Maybe this could be done by extracting the knots positions of the GAM whose knots have been placed automatically by the algorithm, and then use them in the "knots" argument of the GAM in which I need to specify the knots positions.
DUMMY CODE
Note that, in my dummy example, mod_1 (GAM without specifying knots positions) has the highest AIC, while mod_2 (GAM providing knots positions) has the lowest AIC. Please pretende that it is the other way around and the lowest AIC belongs to mod_1 (or edit my code).
library(tidyverse)
library(mgcv)

set.seed(248)

# Produce a table of a dummy, ecological response variable and five predictor variables
# The "soil" variable is a control (categorical) variable, and the x,y coordinates will be
# controlling for spatial variation
eco_data <- tibble(response = round(runif(10000, -0.3, 1.3), 2),
                   soil = as.factor(sample(c("sandy", "clay", "silt"), 10000, T)),
                   rain = sample(800:3000, 10000, T),
                   temp = sample(9:19, 10000, T),
                   xcor = sample(300000:500000, 10000, T),
                   ycor = sample(-450000:-400000, 10000, T))

# GAM without specifying the positions of the knots: AIC = 13052.03
mod_1 <- bam(response ~ soil +
                        s(rain, k = 100, bs = "tp") +
                        s(temp, k = 10, bs = "tp") +
                        s(xcor, ycor, k = 100, bs = "tp"),
             data = eco_data,
             method = "fREML",
             control = gam.control(scalePenalty = FALSE),
             discrete = TRUE)

# Run another model, but specifying the knots
# The positions of the first and last knots will be the min and max values, respectively,
# of the variable, but all other knots should be placed at known locations. This is because
# we want to average the model coefficients with the coefficients of other models,
# so knots must be placed on the same positions in all GAMs
lim_rain <- as.numeric(c(min(eco_data$rain), seq(900, 2900, length = 98), max(eco_data$rain)))
lim_temp <- as.numeric(c(min(eco_data$temp), seq(11, 17, length = 8), max(eco_data$temp)))
lim_xcor <- as.numeric(c(min(eco_data$xcor), seq(320000, 480000, length = 98), max(eco_data$xcor)))
lim_ycor <- as.numeric(c(min(eco_data$ycor), seq(-445000, -405000, length = 98), max(eco_data$ycor)))

# Put all knots into a list
kts <- list(rain = lim_rain, temp = lim_temp, xcor = lim_xcor, ycor = lim_ycor)

# GAM specifying the knots of the four smooth terms: AIC = 12902.49
mod_2 <- bam(response ~ soil +
                        s(rain, k = 100, bs = "tp") +
                        s(temp, k = 10, bs = "tp") +
                        s(xcor, ycor, k = 100, bs = "tp"),
             data = eco_data,
             method = "fREML",
             knots = kts,
             control = gam.control(scalePenalty = FALSE),
             discrete = TRUE)

# The difference in AICs is large: 149.54 (R-squares are also different)
# In my real dataset, by specifying the positions of the knots, I am fitting a worse model
# (diagnostic plots look bad too)

 A: With the thin plate spline, there is a knot at every unique value or unique combination of values of the covariates involved in the smooth. What mgcv then does, because using such a large and rich basis would wasteful,  is to eigen decompose the full basis and then take the k eigenvectors with the k largest eigenvalues as a new basis. This new basis preserves much of the original basis but is of much lower dimensionality.
When you provide k knots, mgcv skips all the above and uses the basis generated assuming those knots only. This basis will likely be less rich than a basis of size k generated using the eigen decomposition, spanning a space of functions that is much smaller than the eigen decomposed one (all else equal).
As the space of functions representable by the basis where you specified exactly k knots is more limited, the model fit is likely to be lower (unless the true function lies in the space spanned by the basis) and this is what is being reflected in the higher AIC when you specify the knots.
If you want to specify the knots, you'll need to generate a rich basis (many knots, as many as the unique data points in the smooth), but if the covariate values change each time the knot locations and hence the basis will change. If you want to fix the basis, you could perhaps find the limits of the covariate over all replicates and then spread n (where n is the number of observations in a replicate) knots over this interval, where n >> k.
