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
x.1
andx.2
for example and doingtest_data[1:30, ]
drops the empty dimension such thattraining_data
isn't even valid input to thedata
argument togam
. Also note you call assigntest_data1
but refer totest_data
. Finally; don't rip components out of models like this, use the extractors like you did forcoef()
:resid(mod.object, type = "response")
andfitted(mod.object)
respectively. Please do run your code in a clean/new/empty session to check that what you posted works... $\endgroup$