¿Las splines están sobreajustando los datos?

Mi problema : recientemente conocí a un estadístico que me informó que las splines solo son útiles para explorar datos y están sujetas a un sobreajuste, por lo que no son útiles en la predicción. Prefería explorar con polinomios simples ... Como soy un gran admirador de las splines, y esto va en contra de mi intuición, estoy interesado en descubrir qué tan válidos son estos argumentos, y si hay un gran grupo de anti-spline- activistas por ahí?

Antecedentes : cuando creo mis modelos, trato de seguir a Frank Harrell, Estrategias de modelado de regresión (1). Argumenta que las splines cúbicas restringidas son una herramienta válida para explorar variables continuas. También argumenta que los polinomios son pobres para modelar ciertas relaciones como umbrales, logarítmicos (2). Para probar la linealidad del modelo, sugiere una prueba ANOVA para la spline:

$H_0: \beta_2 = \beta_3 = … = \beta_{k-1} = 0$

Busqué en Google para sobreajustar con splines, pero no encontré que fuera muy útil (aparte de las advertencias generales sobre no usar demasiados nudos). En este foro parece haber una preferencia por el modelado de splines, Kolassa , Harrell , gung .

Encontré una publicación de blog sobre polinomios, el demonio del sobreajuste que habla sobre la predicción de polinomios. La publicación termina con estos comentarios:

Hasta cierto punto, los ejemplos presentados aquí son trampas: se sabe que la regresión polinómica es altamente no robusta. Mucho mejor en la práctica es usar splines en lugar de polinomios.

Ahora esto me llevó a comprobar cómo funcionarían las splines con el ejemplo:

library(rms)
p4 <- poly(1:100, degree=4)
true4 <- p4 %*% c(1,2,-6,9)
days <- 1:70

set.seed(7987)
noise4 <- true4 + rnorm(100, sd=.5)
reg.n4.4 <- lm(noise4[1:70] ~ poly(days, 4))
reg.n4.4ns <- lm(noise4[1:70] ~ ns(days,4))
dd <- datadist(noise4[1:70], days)
options("datadist" = "dd")
reg.n4.4rcs_ols <- ols(noise4[1:70] ~ rcs(days,5))

plot(1:100, noise4)
nd <- data.frame(days=1:100)
lines(1:100, predict(reg.n4.4, newdata=nd), col="orange", lwd=3)
lines(1:100, predict(reg.n4.4ns, newdata=nd), col="red", lwd=3)
lines(1:100, predict(reg.n4.4rcs_ols, newdata=nd), col="darkblue", lwd=3)

legend("top", fill=c("orange", "red","darkblue"), 
       legend=c("Poly", "Natural splines", "RCS - ols"))

Da la siguiente imagen:

En conclusión, no he encontrado mucho que me convenza de reconsiderar splines, ¿qué me estoy perdiendo?

1. FE Harrell, Estrategias de modelado de regresión: con aplicaciones a modelos lineales, regresión logística y análisis de supervivencia, reimpresión de tapa blanda de tapa dura 1ª ed. 2001. Springer, 2010.
2. FE Harrell, KL Lee y BG Pollock, "Modelos de regresión en estudios clínicos: determinación de las relaciones entre predictores y respuesta", JNCI J Natl Cancer Inst, vol. 80, no. 15, págs. 1198–1202, octubre de 1988.

Actualizar

Los comentarios me hicieron preguntarme qué sucede dentro del intervalo de datos pero con curvas incómodas. En la mayoría de las situaciones, no salgo del límite de datos, como lo indica el ejemplo anterior. No estoy seguro de que esto califique como predicción ...

De todos modos, aquí hay un ejemplo en el que creo una línea más compleja que no se puede traducir a un polinomio. Como la mayoría de las observaciones están en el centro de los datos, traté de simular eso también:

library(rms)
cmplx_line <-  1:200/10
cmplx_line <- cmplx_line + 0.05*(cmplx_line - quantile(cmplx_line, .7))^2
cmplx_line <- cmplx_line - 0.06*(cmplx_line - quantile(cmplx_line, .3))^2
center <- (length(cmplx_line)/4*2):(length(cmplx_line)/4*3)
cmplx_line[center] <- cmplx_line[center] + 
    dnorm(6*(1:length(center)-length(center)/2)/length(center))*10

ds <- data.frame(cmplx_line, x=1:200)

days <- 1:140/2

set.seed(1234)
sample <- round(rnorm(600, mean=100, 60))
sample <- sample[sample <= max(ds$x) & 
                     sample >= min(ds$x)]
sample_ds <- ds[sample, ]

sample_ds$noise4 <- sample_ds$cmplx_line + rnorm(nrow(sample_ds), sd=2)
reg.n4.4 <- lm(noise4 ~ poly(x, 6), data=sample_ds)
dd <- datadist(sample_ds)
options("datadist" = "dd")
reg.n4.4rcs_ols <- ols(noise4 ~ rcs(x, 7), data=sample_ds)
AIC(reg.n4.4)

plot(sample_ds$x, sample_ds$noise4, col="#AAAAAA")
lines(x=ds$x, y=ds$cmplx_line, lwd=3, col="black", lty=4)

nd <- data.frame(x=ds$x)
lines(ds$x, predict(reg.n4.4, newdata=ds), col="orange", lwd=3)
lines(ds$x, predict(reg.n4.4rcs_ols, newdata=ds), col="lightblue", lwd=3)

legend("bottomright", fill=c("black", "orange","lightblue"), 
       legend=c("True line", "Poly", "RCS - ols"), inset=.05)

Esto da la siguiente trama:

Actualización 2

Desde esta publicación, publiqué un artículo que analiza la no linealidad de la edad en un gran conjunto de datos. El suplemento compara diferentes métodos y he escrito una publicación de blog al respecto .

El sobreajuste proviene de permitir una clase de modelos demasiado grande. Esto se vuelve un poco complicado con modelos con parámetros continuos (como splines y polinomios), pero si discretiza los parámetros en una cantidad de valores distintos, verá que al aumentar la cantidad de nudos / coeficientes aumentará la cantidad de modelos disponibles exponencialmente . Para cada conjunto de datos hay una spline y un polinomio que se ajusta con precisión, siempre que permita suficientes coeficientes / nudos. Puede ser que una spline con tres nudos sobrepase más que un polinomio con tres coeficientes, pero eso no es una comparación justa.

Si tiene un número bajo de parámetros y un gran conjunto de datos, puede estar razonablemente seguro de que no está sobreajustando. Si desea probar un mayor número de parámetros, puede intentar la validación cruzada dentro de su conjunto de pruebas para encontrar el mejor número, o puede usar un criterio como Longitud mínima de descripción .

EDITAR : Como se solicitó en los comentarios, un ejemplo de cómo se aplicaría MDL. Primero debe lidiar con el hecho de que sus datos son continuos, por lo que no se pueden representar en un código finito. En aras de la simplicidad, segmentaremos el espacio de datos en cuadros de lado y, en lugar de describir los puntos de datos, describiremos los cuadros en los que se encuentran los datos. Esto significa que perdemos algo de precisión, pero podemos hacer que arbitrariamente pequeño, por lo que no importa mucho.$\epsilon$$\epsilon$

Ahora, la tarea es describir el conjunto de datos de la manera más sucinta posible con la ayuda de algún polinomio. Primero describimos el polinomio. Si es un polinomio de enésimo orden, solo necesitamos almacenar coeficientes (n + 1). Nuevamente, necesitamos discretizar estos valores. Después de eso, debemos almacenar primero el valor en una codificación libre de prefijos (para saber cuándo dejar de leer) y luego eln + $n$$n+1$valores paramétricos. Con esta información, un receptor de nuestro código podría restaurar el polinomio. Luego agregamos el resto de la información requerida para almacenar el conjunto de datos. Para cada punto de datos le damos el valor de x, y luego cuántos cuadros arriba o abajo del punto de datos se encuentra fuera del polinomio. Ambos valores los almacenamos en una codificación libre de prefijos para que los valores cortos requieran pocos bits y no necesitemos delimitadores entre puntos. (Puede acortar el código para los valores de x almacenando solo los incrementos entre valores)

El punto fundamental aquí es la compensación. Si elijo un polinomio de orden a (como f (x) = 3.4), entonces el modelo es muy simple de almacenar, pero para los valores de y, esencialmente estoy almacenando la distancia a la media. Más coeficientes me dan un polinomio de mejor ajuste (y, por lo tanto, códigos más cortos para los valores de y), pero tengo que gastar más bits describiendo el modelo. El modelo que le proporciona el código más corto para sus datos es el que mejor se ajusta al criterio MDL.

(Tenga en cuenta que esto se conoce como 'MDL crudo', y hay algunas mejoras que puede hacer para resolver varios problemas técnicos).

Los estadísticos han estado discutiendo sobre el ajuste polinómico durante años, y en mi experiencia, todo se reduce a esto:

Las splines son básicamente una serie de diferentes ecuaciones unidas, lo que tiende a aumentar la precisión de los valores interpolados a costa de la capacidad de proyectar fuera del rango de datos. Esto está bien si sabe que sus datos son puros y provienen de una fuente constante y si está tratando de describir la probabilidad de la presencia de diferentes valores dentro de su rango de valores. Sin embargo, generalmente no aprendemos tanto sobre los fundamentos teóricos que impulsan nuestros datos, ya que una nueva spline comienza cuando la antigua spline deja de describir con precisión los datos. Esto hace que la predicción de valores fuera de nuestros datos sea casi inútil.

Ahora, las splines no son únicas a este respecto. Las funciones polinomiales en realidad sufren el mismo problema si solo ajustamos los datos y no usamos un marco teórico para elegir las variables. Aquellos que tienen una teoría bien formada que impulsa qué variables permitir que varíen y cuánto confiarán más en la capacidad de una función polinómica compleja para extrapolar predicciones fuera de los datos.

Sin embargo, muchos estadísticos están trabajando con datos sin la ayuda de un marco teórico preestablecido, y esto empuja a algunas personas hacia polinomios simples. Razonan que una función menos flexible que se ajusta a los datos tiene más probabilidades de predecir con precisión los valores fuera de los datos, porque es menos probable que la función se vea influida por anomalías dentro de los datos. Si bien he tenido conversaciones sobre esto con personas que prefieren polinomios simples, nunca tuve la sensación de un grupo anti-spline. Parece que los polinomios simples hacen que algunos estadísticos se sientan más cómodos al evitar el sobreajuste.

Descargo de responsabilidad

Personalmente, no tiendo a usar splines o polinomios simples con la mayoría de mis datos, porque trabajo en un campo con muchos marcos teóricos preestablecidos. Además, generalmente he observado la recopilación de datos y puedo obtener una comprensión decente de lo que estaba impulsando los resultados. En ese caso, estoy construyendo más de un algoritmo lógico y probando la aptitud del algoritmo, en lugar de probar la aptitud de una función polinómica. Puedes agregar este grano de sal a mi respuesta.