Residual bootstrap using a GAM model with two input variables

Question

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.

Thanks

test_data1<-data.frame(y=c(0.20,0.20,0.21,0.21,0.21,0.20,0.19,0.18,0.16,0.10,
                           0.02,-0.02,0.01,0.03,0.07,0.14,0.22,0.13,0.12,
                           0.16,0.17,0.18,0.18,0.17,0.15,0.15,0.13,0.12,
                           0.10,0.08,0.06,0.04,0.03,0.02,0.03,0.05,0.34,
                           0.13,0.11,0.12), 
                       x.1=c(NA,0.20,0.20,0.21,0.21,0.21,0.20,0.19,0.18,0.16,
                             0.10,0.02,-0.02,0.01,0.03,0.07,0.14,0.22,0.13,
                             0.12,0.16,0.17,0.18,0.18,0.17,0.15,0.15,0.13,
                             0.12,0.10,0.08,0.06,0.04,0.03,0.02,0.03,0.05,
                             0.34,0.13,0.11),
                       x.2=c(NA,NA,0.20,0.20,0.21,0.21,0.21,0.20,0.19,0.18,
                             0.16,0.10,0.02,-0.02,0.01,0.03,0.07,0.14,0.22,
                             0.13,0.12,0.16,0.17,0.18,0.18,0.17,0.15,0.15,
                             0.13,0.12,0.10,0.08,0.06,0.04,0.03,0.02,0.03,
                             0.05,0.34,0.13))
training_data<- test_data1[1:30,]
test_gam<- gam(y~ s(x.1, bs="cr")+ s(x.2, bs="cr"), data=training_data)

N=1000
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

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.)