# Bootstrap yielding counterintuitive results (Update: but not anymore)

I have some experience in using bootstrap methods and I'm back to them after a really long hiatus. However, I'm almost certain I'm doing something wrong and, after a lot of time trying to figure out what it is, I'm sure I won't find out by myself.

I'm going to provide a MWE to see if anyone can help me. My simplest attempt is trying to test:

$\begin{cases} H_0:NOx=\displaystyle \sum_{i=0}^5a_iE^i+\varepsilon \text{ for some }a\in\mathbb R^6\\ H_1:NOx=s(E)+\varepsilon \text{ for some smooth }s\text{ (that's not a polynomial of degree$\leq5$) } \end{cases}$

for continuous RVs $NOx$ and $E$ from the ethanol data at the R package lattice. So the models would be:

require(lattice);data(ethanol)
M0<-lm(NOx~E+I(E^2)+I(E^3)+I(E^4)+I(E^5),data=ethanol)
M1<-mgcv::gam(NOx~s(E),data=ethanol)


My test statistic is the relative difference of RSS:

RSS0<-sum(residuals(M0)^2)


and I approximate its null distribution via wild bootstrap (the gold ratio one, Mammen [1993]) the following way:

adj0<-predict(M0)
res0<-residuals(M0)
sigma0<-sd(res0)
n<-nrow(ethanol)

set.seed(1)
B<-1000;Rstar<-rep(NA,B);ethanolstar<-ethanol
veplus<-res0*(1+sqrt(5))/2
veminus<-res0*(1-sqrt(5))/2
for (b in 1:B){
ii<-rbinom(n,1,(5+sqrt(5))/10)
ethanolstar$NOx<-adj0+veminus*ii+veplus*(1-ii) M0star<-lm(NOx~E+I(E^2)+I(E^3)+I(E^4)+I(E^5),data=ethanolstar) M1star<-mgcv::gam(NOx~s(E),data=ethanolstar) RSS1star<-sum(residuals(M1star)^2) RSS0star<-sum(residuals(M0star)^2) Rstar[b]<-(RSS0star-RSS1star)/RSS1star } cat("p-value:",mean(Rstar>R),"\n")  Finally, I get a p-value of 0.01. Similarly, when I use the simple Gaussian bootstrap, that is, defining: ethanolstar$NOx<-rnorm(n,adj0,sigma0)


at each iteration, I get a p-value of 0.001.

Why I suspect these results to be wrong.

A simple plot suggests that a low p-value shouldn't be expected: And, on top of that, I get really low p-values for a variety of examples in which the null model is correct.

So, (why) is my R code wrong?

• I don't know very much about bootstrap methods, but are there problems caused by mgcv doing model selection inside the simulation? Of your 1000 M1Star models, 750 report a lower "effective degrees of freedom" than M1. The 10 simulated Rstar values bigger than R all occur when M1Star edf is bigger than M1's edf. – Jonny Lomond May 29 '18 at 14:22
• First of all, thanks for the feedback! I'm not using mgcv::gam to do model selectionâ€”I just use it to obtain the corresponding residuals. Anyway, I don't think I'm understanding your comment very well... How are those "lower effective degrees of freedom" making the results go wrong? What do EDF have to do with the problem? – jpz May 29 '18 at 15:20
• Why do you believe that the test should not reject based on the plot? It looks to me like the polynomial fit is horrible in the right part of the plot; that could easily lead to rejection. – guy May 29 '18 at 16:14
• @guy exactly, there is a 20% reduction in residual sum of squares using the GAM model. – AdamO May 29 '18 at 16:18
• (+1) @jpz, I ran your code and as far as I can tell, it seems okay. I tried re-running the Gaussian bootstrap on a single dataset generated under $H_0$ (itself from a Gaussian bootstrap) and got a nominal rejection rate as expected. – half-pass May 29 '18 at 16:56

I think the bootstrap is not wrong. What's wrong is the test statistic you're using. The residual sum of squares is calculated using internal validation. The GAM is overfitting the curve. Polynomials approximations obtained from least squares are known to fit relatively badly in the tails, which is actually where a reasonable amount of the ethanol exposure measures are concentrated. Try modifying your approach to use split-sample validation.

I get a much more conservative estimate of the incremental predictive accuracy of the GAM using this suggested modification of your program:

require(lattice)
library(mgcv)

data(ethanol)
set.seed(123)
## 80/20 test/valid
nr <- nrow(ethanol)
split0 <- sample(1:nr, floor(0.8*nr))
split1 <- setdiff(1:nr, split0)

## to generate data under the null
M0 <- lm(NOx~E+I(E^2)+I(E^3)+I(E^4)+I(E^5),data=ethanol[split0,])
res0 <- residuals(M0)

fitntest <- function(ethanol, split0, split1) {
M0 <- lm(NOx~E+I(E^2)+I(E^3)+I(E^4)+I(E^5),data=ethanol[split0,])
M1 <- gam(NOx~s(E),data=ethanol[split0,])
RSS0 <- sum({ethanol[split1, 'NOx']-predict(M0, newdata=ethanol[split1, ])}^2)
RSS1 <- sum({ethanol[split1, 'NOx']-predict(M1, newdata=ethanol[split1, ])}^2)
ethanol$NOx <- predict(M0, newdata=ethanol) + rnorm(nr, 0, sd(res0)) fitntest(ethanol, split0, split1) }) cat("p-value:",mean(Rstar>R),"\n")  • Thank you very much for your answer, @AdamO! I've used your approach with a few examples and it seems to work fine! But I don't know if I'm understanding the formal justification of the need for split-sample validation. Can anybody tell me if what I say below is right? – jpz May 29 '18 at 21:01 • In my OP I use the whole sample to find the best$\hat a\in\mathbb R^6$by looking at the residuals of M0, and the same for$\hat s$in M1. Then I compute the RSS's with the same data, which means that these RSS's will be as low as possible (overfitting). So in the end it's the same as fitting a model with some data and affirming that the model's great because it fits the data very well, which is the typical circular argument that motivates cross-validation. – jpz May 29 '18 at 21:01 • @jpz the reason is that internal validation is too optimistic and so favors the GAM which has a tendency to overfit relative to a 5 degree polynomial. Using$k\$-fold cross validation is yet another way to demote the "optimistm" of the GAM. By all means consider it! – AdamO May 29 '18 at 21:04
• @jpz (Use backticks  to monospace code in comments). I think that's right. Use the split0 and split1` vector to select residuals from the respective train/valid samples. Try it and report back? – AdamO May 30 '18 at 15:41