I've got a GAM model that is in the form of $$Y=f(x_1)+f(x_2)$$ and I would like to perform a residual bootstrap with replacement. Is there any good place where I can see a coded example that does residual bootstrap on a GAM model in R? All the ones I have found so far are done on linear models or generalized linear models.

I've tried creating a function myself to try to create a residual bootstrap but my final plot of the bootstraps produces a high number of frequencies. For example, below I run $N=1000$ bootstraps and in return, I get frequencies larger than the actual number of bootstraps which is not what I was expecting. I've posted the code of my work at the very bottom as a reference.


Example of bootstrap plot.

training_data<- test_data1[1:30,]
test_gam<- gam(y~ s(x.1, bs="cr")+ s(x.2, bs="cr"), data=training_data)

BootstarpFromResiduals<- function(mod.object= test_gam, dat= training_data){
  resids= mod.object$residuals # Extracts residuals from the model
  fittedValues= mod.object$fitted.values
  matr<- model.matrix(mod.object)
  # generating new values for each y[i], by adding bootstrapped resids
  # to fitted values.
  Y<- fittedValues+ sample(resids, length(resids), replace=T)
  # Using model.matrix for the predictors
  model.boot<- gam(Y~ 0+matr, data=dat) # refit model with new Y values
  coef(model.boot) # Extract the coefficients

residual.boot.N<- t(replicate(N, BootstarpFromResiduals()))
# Plots
hist(residual.boot.N) # Plots the bootstrapped residual
  • $\begingroup$ If your observations have no variation around their expected value, as your model indicates, you don't need GAM $\endgroup$
    – Glen_b
    Mar 23, 2022 at 1:07
  • 1
    $\begingroup$ Your code isn't reproducible, there's no x.1 and x.2 for example and doing test_data[1:30, ] drops the empty dimension such that training_data isn't even valid input to the data argument to gam. Also note you call assign test_data1 but refer to test_data. Finally; don't rip components out of models like this, use the extractors like you did for coef(): resid(mod.object, type = "response") and fitted(mod.object) respectively. Please do run your code in a clean/new/empty session to check that what you posted works... $\endgroup$ Mar 23, 2022 at 9:20
  • $\begingroup$ Hello to you both. I'll fix the error with the test data as I hadn't even realized that I had such a mix-up on it. Thanks for the response. $\endgroup$ Mar 23, 2022 at 14:36
  • $\begingroup$ @GavinSimpson thanks for the comment, I cleared the full session and now the code should be completely reproducible. I'll also do the adjustments that you recommend to make the code run more smoothly. $\endgroup$ Mar 23, 2022 at 15:44
  • $\begingroup$ @Warhawk1987 That's fine, but as my answer explains, you generated 19000 coefficient values and drew a histogram of those values explaining the frequencies above 1000 that you were questioning. You'll need to explain what you wanted to do as what you asked R to do, it did, but it doesn't make any sense to ask R to do this. $\endgroup$ Mar 24, 2022 at 10:06

1 Answer 1


You asked for 1000 parametric bootstrap samples and each sample returns 1 + 9 + 9 = 19 coefficients. This is returned from replicate() as a matrix so the dimension of that will be 1000 x 19 or that transposed depending on how replicate() arranges the output.

Regardless, you have a matrix with 19000 elements and you plotted this matrix using hist(), which just treated the matrix as a numeric vector (which it literally is in R) and hence drew a histrogram of 19000 values, leading to the frequencies you are questioning.

(What did you expect hist() to do with a matrix? Or perhaps better, what are you trying to do? Expand on your problem requirements in your question.)


Your Answer

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

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