Regresión polinomial en R (paso a paso)

La regresión polinomial es una técnica que podemos utilizar cuando la relación entre una variable predictora y una variable de respuesta no es lineal.

Este tipo de regresión toma la forma:

Y = β ₀ + β ₁ X + β ₂ X ² +… + β _h X ^h + ε

donde h es el «grado» del polinomio.

Este tutorial proporciona un ejemplo paso a paso de cómo realizar una regresión polinomial en R.

Paso 1: crear los datos

Para este ejemplo, crearemos un conjunto de datos que contiene la cantidad de horas estudiadas y la puntuación del examen final para una clase de 50 estudiantes:

#Haga que este ejemplo sea reproducible
set.seed (1)

#create dataset
 df <- data.frame (hours = runif (50, 5, 15), score = 50)
df $ puntaje = df $ puntaje + df $ horas ^ 3/150 + df $ horas * runif (50, 1, 2)

#ver las primeras seis filas del encabezado de datos
 (datos)

      puntaje de horas
1 7.655087 64.30191
2 8.721239 70.65430
3 10.728534 73.66114
4 14.082078 86.14630
5 7.016819 59.81595
6 13.983897 83.60510

Paso 2: Visualice los datos

Antes de ajustar un modelo de regresión a los datos, primero creemos un diagrama de dispersión para visualizar la relación entre las horas estudiadas y la puntuación del examen:

biblioteca (ggplot2)

ggplot (df, aes (x = horas, y = puntuación)) +
  geom_point ()

Podemos ver que los datos muestran una relación cuadrática, lo que indica que la regresión polinomial podría ajustarse mejor a los datos que la regresión lineal simple.

Paso 3: ajuste los modelos de regresión polinomial

A continuación, ajustaremos cinco modelos de regresión polinomial diferentes con grados h = 1… 5 y usaremos k veces la validación cruzada con k = 10 veces para calcular el MSE de prueba para cada modelo:

# aleatoriamente barajar datos
 df.shuffled <- df [ sample ( nrow (df)),]

# definir el número de pliegues que se utilizarán para la validación cruzada de k pliegues
 K <- 10

#definir el grado de polinomios para ajustar el
 grado <- 5

#create k pliegues de igual tamaño pliegues
 <- cut ( seq (1, nrow (df.shuffled)), breaks = K, labels = FALSE )

# crear objeto para contener MSE de modelos
 mse = matriz (datos = NA, nrow = K, ncol = grado)

#Realice una validación cruzada de K-fold
 para (i en 1: K) {
    
    #define los datos de entrenamiento y prueba
     testIndexes <- que (pliegues == i, arr.ind = TRUE )
    testData <- df.shuffled [testIndexes,]
    trainData <- df.shuffled [-testIndexes,]
    
    #use k-fold cv para evaluar modelos
     para (j en 1: grado) {
        fit.train = lm (puntuación ~ poli (horas, j), datos = trainData)
        fit.test = predecir (fit.train, newdata = testData)
        mse [i, j] = mean ((fit.test-testData $ score) ^ 2)
    }
}

#find MSE para cada grado 
 colMeans (mse)

[1] 9.802397 8.748666 9.601865 10.592569 13.545547

Desde el resultado podemos ver el MSE de prueba para cada modelo:

• Prueba MSE con grado h = 1: 9.80
• Prueba MSE con grado h = 2: 8,75
• Prueba MSE con grado h = 3: 9.60
• Prueba MSE con grado h = 4: 10.59
• Prueba MSE con grado h = 5: 13.55

El modelo con el MSE de prueba más bajo resultó ser el modelo de regresión polinomial con grado h = 2.

Esto coincide con nuestra intuición del diagrama de dispersión original: un modelo de regresión cuadrática se ajusta mejor a los datos.

Paso 4: analizar el modelo final

Por último, podemos obtener los coeficientes del modelo de mejor rendimiento:

# se ajusta al mejor modelo 
mejor = lm (puntuación ~ poli (horas, 2, sin procesar = T ), datos = gl)

#ver resumen del mejor modelo
resumen (mejor)

Llamada:
lm (fórmula = puntuación ~ poli (horas, 2, sin procesar = T), datos = gl)

Derechos residuales de autor:
    Mín. 1T Mediana 3T Máx. 
-5,6589 -2,0770 -0,4599 2,5923 4,5122 

Coeficientes:
                         Estimar Std. Valor t de error Pr (> | t |)    
(Intercepción) 54.00526 5.52855 9.768 6.78e-13 ***
poli (horas, 2, crudo = T) 1 -0.07904 1.15413 -0.068 0.94569    
poli (horas, 2, crudo = T) 2 0.18596 0.05724 3.249 0.00214 ** 
---
Signif. códigos: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0,1 pulg. 1

De la salida podemos ver que el modelo ajustado final es:

Puntuación = 54.00526 – .07904 * (horas) + .18596 * (horas) ²

Podemos usar esta ecuación para estimar la puntuación que recibirá un estudiante en función de la cantidad de horas que estudió.

Por ejemplo, se espera que un estudiante que estudia durante 10 horas reciba una puntuación de 71,81 :

Puntaje = 54.00526 – .07904 * (10) + .18596 * (10) ² = 71.81

También podemos trazar el modelo ajustado para ver qué tan bien se ajusta a los datos sin procesar:

ggplot (df, aes (x = horas, y = puntuación)) +
          geom_point () +
          stat_smooth (método = ' lm ', fórmula = y ~ poli (x, 2), tamaño = 1) +
          xlab (' Horas estudiadas ') +
          ylab (' Puntuación ')

Puede encontrar el código R completo utilizado en este ejemplo aquí .

• https://r-project.org
• https://www.python.org/
• https://www.stata.com/