Bootstrap dando resultados poco intuitivos (Actualización: pero ya no)

Question

Tengo algo de experiencia en el uso de métodos bootstrap y yo estoy de vuelta, después de una muy larga ausencia. Sin embargo, estoy casi segura de que estoy haciendo algo mal y, después de un montón de tiempo tratando de averiguar lo que es, estoy seguro de que no voy a averiguar por mí mismo.

---

Me voy a dar una MWE a ver si alguien me puede ayudar. Mi más simple intento es tratar de probar:

$\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 %#%#%) } \end{casos}$\leq5$NOx$

for continuous RVs $E$ a partir del etanol de datos en el paquete de R celosía. Por lo que los modelos sería:

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)

---

Mi estadístico de prueba es la diferencia relativa de RSS:

RSS0<-sum(residuals(M0)^2)
RSS1<-sum(residuals(M1)^2)
R<-(RSS0-RSS1)/RSS1

y yo aproximado de su nula distribución a través de wild bootstrap (el oro relación de uno, Mammen [1993]) de la siguiente manera:

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

---

Finalmente, obtener un p-valor de 0.01.

Del mismo modo, cuando uso el simple Gaussiano bootstrap, que es, la definición de:

ethanolstar$NOx<-rnorm(n,adj0,sigma0)

en cada iteración, puedo obtener un p-valor de 0,001.

---

Por qué tengo la sospecha de que estos resultados sean mal.

Una trama muy simple sugiere que un bajo p-valor no debe ser esperado: Y, encima de eso, me pongo muy bajos de p-valores para una variedad de ejemplos en los que el modelo nulo es correcta.

---

Así, (por qué) es mi código R mal?

Answer 1

Creo que el bootstrap no es malo. Lo que está mal es la estadística de prueba de que usted está utilizando. La suma de cuadrados residual se calcula mediante la validación interna. El GAM es el sobreajuste de la curva. Los polinomios de aproximaciones obtenidas a partir de los mínimos cuadrados son conocidos para que se ajuste relativamente mal en las colas, que es donde en realidad una cantidad razonable de la exposición a etanol medidas se concentran. Trate de modificar su estrategia de división de la muestra de validación.

Puedo obtener una estimación conservadora de la incrementales de la exactitud de predicción de la GAM el uso de este sugirió la modificación de su programa:

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