Bootstrap yielding counterintuitive results (Update: but not anymore)

Question

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)
RSS1<-sum(residuals(M1)^2)
R<-(RSS0-RSS1)/RSS1

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?

Answer 1

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)
  (RSS0-RSS1)/RSS1
}

R <- fitntest(ethanol, split0, split1)
Rstar <- replicate(1000, {
  ethanol$NOx <- predict(M0, newdata=ethanol) + rnorm(nr, 0, sd(res0))
  fitntest(ethanol, split0, split1)
})
cat("p-value:",mean(Rstar>R),"\n")