3
$\begingroup$

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)
$\endgroup$
3
  • 1
    $\begingroup$ I cannot reproduce these values you are reporting. In particular I get AIC(mod_1) and AIC(mod_2) as 12945.48 and 12944.89 respectively. I am using R version 4.0.2 and mgcv 1.8-33. One of the reasons for the AIC difference observed (here) is that we have different degrees of freedom for the two models; mod_2 seems to have somewhat "fewer" (12.70 instead of 13.31). Notice that the DF in GAM need to account for the smoothing and the uncertainty around it. (Interesting question: +1) $\endgroup$
    – usεr11852
    Sep 5, 2020 at 12:44
  • $\begingroup$ I've updated my mgcv package to version 1.8-33 (it was 1.8.-31), but my R version is the same as yours. Now, I'm getting the same results as you. However, the problem still remains in the real data, where the model without specifying the positions of the knots has a far better AIC than the model in which I am providing the positions of knots (delta AIC = 1442.47). The DF are also different (mod_1 = 213.11 vs. mod_2 = 197.49). $\endgroup$
    – Darius
    Sep 6, 2020 at 10:50
  • $\begingroup$ I guess the difference in delta AIC is the positions of knots: in mod_1, the algorithm chooses the best positions, whereas in mod_2, I'm choosing equal-distance positions. My question of whether it would be possibe to achieve a similar model fit (delta AIC < 2) between the two models is still bothering me. $\endgroup$
    – Darius
    Sep 6, 2020 at 10:52

1 Answer 1

3
$\begingroup$

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.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.