Ajuste de Regresión Polinómica en R

. (Puede informar de un problema sobre el contenido en esta página aquí) ¿Desea compartir su contenido en R-bloggers? haga clic aquí si tiene un blog, o aquí si no.

Una relación lineal entre dos variables x y y es una de las suposiciones más comunes, efectivas y fáciles de hacer al tratar de averiguar su relación. A veces, sin embargo, la verdadera relación subyacente es más compleja que eso, y aquí es cuando la regresión polinómica entra en ayuda.

Veamos un ejemplo de economía: Supongamos que desea comprar una cierta cantidad q de un determinado producto. Si el precio unitario es p, entonces usted pagaría una cantidad total y. Este es un ejemplo típico de una relación lineal. El precio total y la cantidad son directamente proporcionales. Para trazarlo escribiríamos algo como esto:

p <- 0.5q <- seq(0,100,1)y <- p*qplot(q,y,type='l',col='red',main='Linear relationship')

La trama se verá así:
linear-relationship

Ahora, esta es una buena aproximación de la verdadera relación entre y y q, sin embargo, al comprar y vender, es posible que deseemos considerar alguna otra información relevante, como: Comprar cantidades significativas es probable que podamos pedir y obtener un descuento, o comprar más y más de un cierto bien que podríamos estar empujando el precio hacia arriba.
Esto puede llevar a un escenario como este en el que el costo total ya no es una función lineal de la cantidad:

y <- 450 + p*(q-10)^3plot(q,y,type='l',col='navy',main='Nonlinear relationship',lwd=3)

Rplot02
Con la regresión polinómica podemos ajustar modelos de orden n > 1 a los datos e intentar modelar relaciones no lineales.

Cómo ajustar una regresión polinómica

En primer lugar, recuerde siempre usar set.seed(n) al generar números pseudo aleatorios. Al hacer esto, el generador de números aleatorios genera siempre los mismos números.

set.seed(20)

Predictor (q). Utilice seq para generar secuencias igualmente espaciadas rápido

q <- seq(from=0, to=20, by=0.1)

Valor para predecir (y):

y <- 500 + 0.4 * (q-10)^3

Se genera algo de ruido y se agrega a la señal real (y):

noise <- rnorm(length(q), mean=10, sd=80)noisy.y <- y + noise

Gráfico de la señal ruidosa:

plot(q,noisy.y,col='deepskyblue4',xlab='q',main='Observed data')lines(q,y,col='firebrick1',lwd=3)

Esta es la gráfica de nuestros datos observados simulados. Los puntos de datos simulados son los puntos azules, mientras que la línea roja es la señal (señal es un término técnico que a menudo se usa para indicar la tendencia general que estamos interesados en detectar).
Rplot03
Nuestro modelo debe ser algo como esto: y = a*q + b*q2 + c*q3 + costo

Ajustémoslo usando R. Al ajustar polinomios, puede usar

model <- lm(noisy.y ~ poly(q,3))

O

model <- lm(noisy.y ~ x + I(X^2) + I(X^3))

Sin embargo, tenga en cuenta que q, I(q^2) y I(q^3) estará correlacionado y las variables correlacionadas pueden causar problemas. El uso de poly() le permite evitar esto produciendo polinomios ortogonales, por lo tanto, voy a usar la primera opción.

summary(model)Call:lm(formula = noisy.y ~ poly(q, 3))Residuals: Min 1Q Median 3Q Max -212.326 -51.186 4.276 61.485 165.960 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 513.615 5.602 91.69 <2e-16 ***poly(q, 3)1 2075.899 79.422 26.14 <2e-16 ***poly(q, 3)2 -108.004 79.422 -1.36 0.175 poly(q, 3)3 864.025 79.422 10.88 <2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 79.42 on 197 degrees of freedomMultiple R-squared: 0.8031,Adjusted R-squared: 0.8001 F-statistic: 267.8 on 3 and 197 DF, p-value: 0 

Utilizando la función confint() podemos obtener los intervalos de confianza de los parámetros de nuestro modelo.
Intervalos de confianza para los parámetros del modelo:

confint(model, level=0.95) 2.5 % 97.5 %(Intercept) 502.5676 524.66261poly(q, 3)1 1919.2739 2232.52494poly(q, 3)2 -264.6292 48.62188poly(q, 3)3 707.3999 1020.65097

Parcela de mobiliario vs residuos. No se debe mostrar un patrón claro en la gráfica de residuos si el modelo se ajusta bien

plot(fitted(model),residuals(model))

Rplot04
En general, el modelo parece un buen ajuste, ya que la R cuadrada de 0.8 indica. Los coeficientes de los términos de primer y tercer orden son estadísticamente significativos como esperábamos. Ahora podemos usar la función predict() para obtener los valores ajustados y los intervalos de confianza con el fin de trazar todo contra nuestros datos.

Valores previstos e intervalos de confianza:

predicted.intervals <- predict(model,data.frame(x=q),interval='confidence', level=0.99)

Añadir líneas a la gráfica existente:

lines(q,predicted.intervals,col='green',lwd=3)lines(q,predicted.intervals,col='black',lwd=1)lines(q,predicted.intervals,col='black',lwd=1)

Añadir una leyenda:

legend("bottomright",c("Observ.","Signal","Predicted"), col=c("deepskyblue4","red","green"), lwd=3)

Aquí está la trama:
Rplot05
Podemos ver que nuestro modelo hizo un trabajo decente al ajustar los datos y, por lo tanto, podemos estar satisfechos con él.

Una advertencia: Los polinomios son herramientas poderosas, pero podrían ser contraproducentes: en este caso, sabíamos que la señal original se generaba utilizando un polinomio de tercer grado, sin embargo, al analizar datos reales, generalmente sabemos poco al respecto y, por lo tanto, debemos ser cautelosos porque el uso de polinomios de alto orden (n > 4) puede conducir a un ajuste excesivo. El ajuste excesivo ocurre cuando su modelo está captando el ruido en lugar de la señal: a pesar de que su modelo está mejorando cada vez más en el ajuste de los datos existentes, esto puede ser malo cuando está tratando de predecir nuevos datos y conducir a resultados engañosos.

Un resumen con el código completo de este ejemplo se puede encontrar aquí.

Gracias por leer este post, deja un comentario a continuación si tienes alguna pregunta.

Deja una respuesta

Tu dirección de correo electrónico no será publicada.