¿Cómo encuentro un valor p de regresión suave de spline / loess?

10

Tengo algunas variables y estoy interesado en encontrar relaciones no lineales entre ellas. Así que decidí ajustar algunas spline o loess, e imprimir bonitos gráficos (ver el código a continuación). Pero también quiero tener algunas estadísticas que me den una idea de la probabilidad de que la relación sea una cuestión de aleatoriedad ... es decir, necesito un valor p general, como el que tengo para la regresión lineal, por ejemplo. En otras palabras, necesito saber si la curva ajustada tiene algún sentido, ya que mi código ajustará una curva a cualquier dato.

x <- rnorm(1000)
y <- sin(x) + rnorm(1000, 0, 0.5)

cor.test(x,y)
plot(x, y, xlab = xlab, ylab = ylab)
spl1 <- smooth.spline(x, y, tol = 1e-6, df = 8)
lines(spl1, col = "green", lwd = 2)

spl2 <- loess(y ~ x)
x.pr <- seq(min(x), max(x), length.out = 100)
lines(x.pr, predict(spl2, x.pr), col = "blue", lwd = 2)
Curioso
fuente

Respuestas:

8

La biblioteca de splines tiene funciones bsy nseso creará una base de splines para usar con la lmfunción, luego puede ajustar un modelo lineal y un modelo que incluye splines y usar la anovafunción para hacer la prueba de modelo completo y reducido para ver si el modelo de spline se ajusta significativamente mejor que el modelo lineal.

Aquí hay un código de ejemplo:

x <- rnorm(1000)
y <- sin(x) + rnorm(1000, 0, 0.5)

library(splines)

fit1 <- lm(y~x)
fit0 <- lm(y~1)
fit2 <- lm(y~bs(x,5))

anova(fit1,fit2)
anova(fit0,fit2)

plot(x,y, pch='.')
abline(fit1, col='red')
xx <- seq(min(x),max(x), length.out=250)
yy <- predict(fit2, data.frame(x=xx))
lines(xx,yy, col='blue')

También puede usar la polyfunción para hacer un ajuste polinómico y probar los términos no lineales como prueba de curvatura.

Para el ajuste loess es un poco más complicado. Hay algunas estimaciones de grados equivalentes de libertad para el parámetro de suavizado de loess que podrían usarse junto con los valores para los modelos lineales y loess para construir y la prueba F. Creo que los métodos basados ​​en pruebas de arranque y permutación pueden ser más intuitivos.R2

Existen técnicas para calcular y trazar un intervalo de confianza para un ajuste de loess (creo que puede haber una forma incorporada en el paquete ggplot2), puede trazar la banda de confianza y ver si una línea recta encajaría dentro de la banda (esto no es un valor p, pero aún así da un sí / no.

Podría ajustar un modelo lineal y tomar los residuos y ajustar un modelo de loess a los residuos como respuesta (y la variable de interés como predictor), si el modelo verdadero es lineal, entonces este ajuste debe estar cerca de una línea plana y reordenar los puntos en relación con el predictor no debería hacer ninguna diferencia. Puede usar esto para crear una prueba de permutación. Ajuste el loess, encuentre el valor predicho más alejado de 0, ahora permute aleatoriamente los puntos y ajuste un nuevo loess y encuentre el punto predicho más alejado de 0, repita varias veces, el valor p es la proporción de valores permutados que están más lejos desde 0 que el valor original.

También puede considerar la validación cruzada como un método para elegir el ancho de banda de loess. Esto no da un valor p, pero un ancho de banda infinito corresponde a un modelo lineal perfecto, si la validación cruzada sugiere un ancho de banda muy grande, entonces eso sugiere que un modelo lineal puede ser razonable, si los anchos de banda más altos son claramente inferiores a algunos de los anchos de banda más pequeños entonces esto sugiere una curvatura definida y lineal no es suficiente.

Greg Snow
fuente
Gracias Greg! Creo que el primer párrafo suena como el camino a seguir, excepto que no estoy interesado en la comparación con el modelo lineal, solo para ver si la spline lo explica o no. ¿Podría proporcionar algún código o más punteros concretos sobre cómo probar la spline con anova? He estado mirando las funciones bs y ns, pero no soy tan bueno en estadística como para poder inventarlo yo mismo.
Curioso
R2R2
anovaR2R21-R2R21-R2
Greg, gracias! 1) ¿Podría explicar qué está lm(y~bs(x,5))haciendo y por qué no lm(y~I(bs(x,5)))? Estoy bastante confundido con esta llamada porque el resultado de bs (x, 5) no es una variable ... 2) ¿Entiendo correctamente que el valor p que estoy buscando es el resultado anova(fit0,fit2)?
Curioso
1
XX2X3bsXlm