Dificultad para probar linealidad en regresión

En Modelización estadística: las dos culturas, Leo Breiman escribe

La práctica actual aplicada es verificar el ajuste del modelo de datos utilizando pruebas de bondad de ajuste y análisis residual. En un momento, hace algunos años, configuré un problema de regresión simulada en siete dimensiones con una cantidad controlada de no linealidad. Las pruebas estándar de bondad de ajuste no rechazaron la linealidad hasta que la no linealidad fuera extrema.

Breiman no da los detalles de su simulación. Hace referencia a un artículo que, según él, da una justificación teórica para su observación, pero el documento no está publicado.

¿Alguien ha visto un resultado de simulación publicado o un documento teórico para respaldar la afirmación de Brieman?

Creé una simulación que respondería a la descripción de Breiman y encontré solo lo obvio: el resultado depende del contexto y de lo que se entiende por "extremo".

Se podría decir mucho, pero permítanme limitarlo a un solo ejemplo realizado mediante un Rcódigo fácilmente modificable para que los lectores interesados ​​lo utilicen en sus propias investigaciones. Este código comienza configurando una matriz de diseño que consta de valores independientes distribuidos aproximadamente de manera uniforme que son aproximadamente ortogonales (para que no entremos en problemas de multicolinealidad). Calcula una única interacción cuadrática (es decir, no lineal) entre las dos primeras variables: este es solo uno de los muchos tipos de "no linealidades" que podrían estudiarse, pero al menos es común y bien entendido. Luego estandariza todo para que los coeficientes sean comparables:

set.seed(41)
p <- 7                                            # Dimensions
n <- 2^p                                          # Observations
x <- as.matrix(do.call(expand.grid, lapply(as.list(1:p), function(i) c(-1,1))))
x <- x + runif(n*p, min=-1, max=1)
x <- cbind(x, x.12 = x[,1]*x[,2])                 # The nonlinear part
x <- apply(x, 2, function(y) (y - mean(y))/sd(y)) # Standardization

Para el modelo OLS base (sin no linealidad) debemos especificar algunos coeficientes y la desviación estándar del error residual. Aquí hay un conjunto de coeficientes unitarios y un SD comparable:

beta <- rep(c(1,-1), p)[1:p]
sd <- 1

Para ilustrar la situación, aquí hay una iteración codificada de la simulación. Genera la variable dependiente, resume sus valores, muestra la matriz de correlación completa de todas las variables (incluida la interacción) y muestra una matriz de diagrama de dispersión. Luego realiza la regresión OLS. A continuación, el coeficiente de interacción de es sustancialmente más pequeño que cualquiera de los otros coeficientes (todos iguales a o ), por lo que sería difícil llamarlo "extremo":$1/4$$1$$-1$

gamma = 1/4          # The standardized interaction term
df <- data.frame(x)
df$y <- x %*% c(beta, gamma) + rnorm(n, sd=sd)
summary(df)
cor(df)*100
plot(df, lower.panel=function(x,y) lines(lowess(y~x)), 
     upper.panel=function(x,y) points(x,y, pch=".", cex=4))
summary(lm(df$y ~ x))

En lugar de pasar por toda la salida aquí, echemos un vistazo a estos datos utilizando la salida del plotcomando:

Las trazas de lowess en el triángulo inferior muestran esencialmente ninguna relación lineal entre la interacción ( x.12) y la variable dependiente ( y) y relaciones lineales modestas entre las otras variables y y. Los resultados de OLS lo confirman; la interacción es escasamente significativa:

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.0263     0.0828    0.32    0.751    
xVar1         0.9947     0.0833   11.94   <2e-16 ***
xVar2        -0.8713     0.0842  -10.35   <2e-16 ***
xVar3         1.0709     0.0836   12.81   <2e-16 ***
xVar4        -1.0007     0.0840  -11.92   <2e-16 ***
xVar5         1.0233     0.0836   12.24   <2e-16 ***
xVar6        -0.9514     0.0835  -11.40   <2e-16 ***
xVar7         1.0482     0.0835   12.56   <2e-16 ***
xx.12         0.1902     0.0836    2.27    0.025 *  

Tomaré el valor p del término de interacción como una prueba de no linealidad: cuando este valor p sea suficientemente bajo (puede elegir qué tan bajo), habremos detectado la no linealidad.

(Aquí hay una sutileza acerca de lo que estamos buscando exactamente. En la práctica, podríamos tener que examinar todas las interacciones cuadráticas posibles 7 * 6/2 = 21, así como quizás 7 términos cuadráticos más, en lugar de centrarnos en un solo término como se hace aquí. Nos gustaría hacer una corrección para estas 28 pruebas interrelacionadas. No hago esta corrección explícitamente aquí, porque en su lugar visualizo la distribución simulada de los valores p. Puede leer las tasas de detección directamente desde los histogramas al final según sus umbrales de importancia).

Pero no hagamos este análisis solo una vez; hagámoslo muchas veces, generando nuevos valores de yen cada iteración de acuerdo con el mismo modelo y la misma matriz de diseño. Para lograr esto, utilizamos una función para llevar a cabo una iteración y devolver el valor p del término de interacción:

test <- function(gamma, sd=1) {
  y <- x %*% c(beta, gamma) + rnorm(n, sd=sd)
  fit <- summary(lm(y ~ x))
  m <- coef(fit)
  n <- dim(m)[1]
  m[n, 4]
}

Elijo presentar los resultados de la simulación como histogramas de los valores p, variando el coeficiente estandarizado gammadel término de interacción. Primero, los histogramas:

h <- function(g, n.trials=1000) {
  hist(replicate(n.trials, test(g, sd)), xlim=c(0,1), 
       main=toString(g), xlab="x1:x2 p-value")
}
par(mfrow=c(2,2)) # Draw a 2 by 2 panel of results

Ahora para hacer el trabajo. Se necesitan unos segundos para 1000 ensayos por simulación (y cuatro simulaciones independientes, comenzando con el valor dado del término de interacción y reduciéndolo a la mitad cada vez):

temp <- sapply(2^(-3:0) * gamma, h)

Los resultados:

Leyendo hacia atrás desde la parte inferior derecha, estos gráficos muestran que para esta matriz de diseño x, para esta desviación estándar de errores sdy para estos coeficientes estandarizados beta, OLS puede detectar una interacción estandarizada de (solo un cuarto del tamaño de los otros coeficientes) ) de manera confiable, más del 80% del tiempo (utilizando un umbral del 5% para el valor p - recuerde la breve discusión sobre la corrección de comparaciones múltiples, que ahora estoy ignorando); a menudo puede detectar un tamaño de interacción de (aproximadamente el 20% del tiempo); a veces detectará una interacción de tamaño$1/4$$1/8$$1/16$, y realmente no puede identificar interacciones más pequeñas. Aquí no se muestra un histograma gammaigual a , que muestra que incluso cuando se corrigen comparaciones múltiples, es casi seguro que se detecta una interacción cuadrática de este tamaño.$1/2$

$1/32$$1/4$xsdbetasd

En resumen, una simulación como esta puede probar lo que quiera si simplemente lo configura e interpreta de la manera correcta. Eso sugiere que el estadístico individual debe realizar sus propias exploraciones, adecuadas a los problemas particulares que enfrentan, a fin de llegar a una comprensión personal y profunda de las capacidades y debilidades de los procedimientos que están utilizando.

No estoy seguro de que dé una respuesta final a la pregunta, pero le daría un vistazo a esto . Especialmente el punto 2. Véase también la discusión en el apéndice A2 del documento .