1 INTRODUCCIÓN: ¿QUE ES LA REGRESIÓN LINEAL?

Es un procedimiento estadístico para crear un modelo matemático (una fórmula) que relacione los valores de una variable \(y\) (variable dependiente) en función de los valores de una variable \(x\) (variable independiente). Es el más sencillo y popular para relacionar y predecir variables cuantitativas continuas.

Este modelo puede utilizarse para predecir el futuro comportamiento de la variable dependiente, en caso de disponer de nuevos valores de la variable independiente.

El modelo matemático puede expresarse con la fórmula \(y = ax + b + e\), en el que:

Figura 1: Cálculo de los términos \(a\) y \(b\) según el método de mínimos cuadrados

Los término \(a\) y \(b\) pueden calcularse según diferentes métodos. El más simple es el de mínimos cuadrados (aunque existan otros como el de máxima verosimilitud), que consiste en encontrar aquellos valores de \(a\) y \(b\) que hagan mínima la suma de los cuadrados de las desviaciones de las observaciones respecto de la recta que representa el modelo, ya que en este caso, la relación entre ambas variables es representada por una línea recta.

2 OBTENCIÓN DE UN MODELO BÁSICO DE REGRESIÓN LINEAL

Supongamos las siguientes variables:

x_ejemplo1 <- c(0.8,1,1.2,1.3)
y_ejemplo1 <-c(1,2,3,5)

Figura 2: Ejemplo de regresión lineal entre dos variables desarrollado

media_x_ejemplo1 <- mean(x_ejemplo1)
media_y_ejemplo1 <- mean(y_ejemplo1)

n_ejemplo1 <- length(x_ejemplo1)                                # Número de casos
x.y_ejemplo1 <- x_ejemplo1 * y_ejemplo1                         # Columna producto X*Y
suma.x.y_ejemplo1 <- sum(x.y_ejemplo1)                          # Suma columna X*Y
cov.x.y_ejemplo1 <- sum(x.y_ejemplo1)/n_ejemplo1 - media_x_ejemplo1 * media_y_ejemplo1
x_cuadrado_ejemplo1 <- x_ejemplo1^2                             # Columna con los cuadrados de X
suma.x_cuadrado_ejemplo1 <- sum(x_cuadrado_ejemplo1)

varianza_x_ejemplo1 <- (suma.x_cuadrado_ejemplo1 / n_ejemplo1) - media_x_ejemplo1^2
# Parámetro A
A_ejemplo1 <- cov.x.y_ejemplo1 / varianza_x_ejemplo1
A_ejemplo1

# Parámetro B
B_ejemplo1 <- media_y_ejemplo1 - A_ejemplo1 * media_x_ejemplo1
B_ejemplo1

Por lo tanto, nuestro modelo de regresión tomaría la forma de \(y = 7.2x - 5.084746\).

rm(list=ls())

3 REGRESIÓN LINEAL EN R

La función básica para obtener un modelo de regresión lineal es lm(). Su sintaxis es la siguiente:

lm(formula, data, …)

donde:

Tomando como referencia el ejemplo elaborado manualmente la sintaxis sería la siguiente:

x_ejemplo1 <- c(0.8,1,1.2,1.3)
y_ejemplo1 <-c(1,2,3,5)
lm(y_ejemplo1 ~ x_ejemplo1)

En la salida del procedimiento, los coefficents equivalen a los parámetros a y b.

4 ELABORACIÓN DE UN MODELO COMPLETO DE REGRESIÓN LINEAL EN R.

Aunque parezca sencillo, la elaboración de un modelo de regresión requiere superar varios pasos:

En este ejemplo recurriremos al análisis de algunas variables del conjunto de datos mtcars. En principio, nos interesa saber si existe relación entre dos variables:

Podemos asumir que el consumo de un vehículo (variable dependiente) vendrá parcialmente determinado por su peso (variable independiente) de tal manera que a más peso recorrerá menos porque consumirá más. Por lo tanto, mpg será \(y\) y wt será \(x\).

data("mtcars")                    # Cargamos el fichero directamente del conjunto de ficheros disponibles en R
str(mtcars)

4.1 Supuestos del modelo de regresión lineal

Para poder utilizar un modelo de regresión lineal es necesario que se cumplan con los siguientes supuestos:

  • Linealidad: que la relación entre las variables sea lineal.

  • Normalidad: que las variables tengan una distribución simétrica y gaussiana.

La observación de la nube de puntos de un gráfico de dispersión puede ayudar a confirmar si existe una relación lineal entre ambas variables,

plot(mtcars$mpg ~ mtcars$wt)                         # dependiente ~ independiente
plot(mtcars$wt, mtcars$mpg)                          # independiente, dependiente

La condición de normalidad de ambas variables puede ser evaluada de diversas maneras

  • Visualmente mediante histogramas. A continuación, se representan los histogramas de ambas variables, sobre los que se superpondrá una curva normal, calculada a partir del valor medio y de la desviación típica de cada variable
media_mpg <- mean(mtcars$mpg)                             # Cálculo de la media
desv_tipica_mpg <- sd(mtcars$mpg)                         # Cálculo de la desviación típica
hist(mtcars$mpg, prob=TRUE,                               # Histograma
     xlab="mpg", ylim=c(0, 0.10), 
     main="Curva normal sobre el histograma")
curve(dnorm(x, mean=media_mpg, sd=desv_tipica_mpg),       # Curva normal
      col="darkblue", lwd=2, add=TRUE)
media_wt <- mean(mtcars$wt)                                   # Cálculo de la media
desv_tipica_wt <- sd(mtcars$wt)                         # Cálculo de la desviación típica
hist(mtcars$wt, prob=TRUE,                                    # Histograma
     xlab="wt", ylim=c(0, 0.60), 
     main="Curva Normal")
curve(dnorm(x, mean=media_wt, sd=desv_tipica_wt),       # Curva normal
      col="darkblue", lwd=2, add=TRUE)
  • Gráfico Q-Q. Es un diagrama de puntos que dibuja las funciones de distribución acumuladas. En caso de que provengan de la misma distribución (en este caso la distribución normal), los puntos debería aparecer alineados. Para facilitar esta visualización, se puede superponer una línea recta.
qqnorm(mtcars$mpg,                                        # Gráfico Q-Q
       pch = 1, 
       frame = FALSE)               
qqline(mtcars$mpg,                                        # Línea recta de referencia añadida
       col = "steelblue",                                 # Color de la línea
       lwd = 2)                                           # Anchura de la línea

El paquete car incorpora la función qqPlot que elabora directamente este tipo de gráfico. Este gráfico incorpora las bandas de confianza, que permiten detectar qué casos pueden ser considerados atípicos.

# install.packages("car")
library("car")
qqPlot(mtcars$mpg)              # Dos puntos aparecen fuera del intervalo de confianza (la zona entre las lineas a rayas):

  • Otra posibilidad para evaluar esta condición es aplicar el test de Shapiro-Wilks. Recordemos que:

    • Hipótesis H0: cuando p-value es mayor que \(\alpha\) = 0,05 -> aceptar H0 -> la variable presenta una distribución normal

    • Hipótesis H1: cuando p-value es menor que \(\alpha\) = 0,05 -> rechazar H0 -> la variable presenta una distribución no normal.

shapiro.test(mtcars$mpg)           
shapiro.test(mtcars$wt)

De acuerdo con lo expuesto, ambas variables siguen aproximadamente una distribución normal, por lo que, en teoría, podríamos someter ambas variables a un análisis de regresión.

rm(desv_tipica_mpg, desviacion_tipica_wt, media_mpg, media_wt)

4.2 Obtención del modelo de regresión.

La función lm() crea un objeto de R que se denominará modelo_regresion_1. Éste es una lista conteniendo la información relevante sobre el análisis realizado.

modelo_regresion_1 <- lm(mpg ~ wt,                          # Variable dependiente ~ Variable independiente
                         data = mtcars)                                         
print(modelo_regresion_1)

De acuerdo con los resultados, la ecuación de la recta de mínimos cuadrados es: \(y= -0.5.344x + 37.285\).

La información contenida en el objeto modelo_regresion_1 es muy abundante. Se puede utilizar la función attributes() para mostrar esa información.

attributes(modelo_regresion_1)

Es posible acceder a cada uno de los elementos individuales de esa lista usando el operador $ entre el objeto y el elemento de la lista. En este caso, los parámetros \(a\) y \(b\).

modelo_regresion_1$coefficients                     # Parámetros $a$ y $b$ de la ecuación

También pueden visualizarse los valores ajustados (valores teóricos de mgp al aplicar el modelo, también llamados predicciones).

modelo_regresion_1$fitted.values

También pueden visualizarse los resíduos, es decir, la diferencia entre el valor real de \(y\) y su valor teórico, obtenido a partir del modelo de regresión. Cuanto más pequeños sean mejor será el ajuste del modelo a los datos y más acertadas serán las predicciones que se realicen a partir de dicho modelo.

modelo_regresion_1$residuals

Los valores ajustados y los residuales también se pueden recuperar usando las funciones fitted() y residuals() aplicadas al objeto resultante de la función lm().

fitted(modelo_regresion_1)
residuals(modelo_regresion_1)

4.3 Resumen del modelo.

El resumen del modelo de regresión se obtiene aplicando la función summary.

summary(modelo_regresion_1)

Los resultados de esta función son los siguientes:

  • Call: muestra la función usada para calcular el modelo de regresión.

  • Residuals: ofrece una rápida visión de la distribución de los resíduos (diferencia entre el valor real de \(y\) y el calculado), cuya media es, por definición, 0. Un buen modelo es aquel en el que la diferencia entre la media y la mediana es próxima a 0, y los valores máximos y mínimos tienen valores absolutos similares.

  • Coefficients: valor de la pendiente y origen del modelo y su significación estadística. Por lo tanto, en este ejemplo la ecuación de la recta de mínimos cuadrados es: \(y = -5.34x + 37.28\).

  • Residuals Standard Error (RSE): métrica que informa de cómo se ajusta el modelo a los datos originales. Cuanto más bajo sea el valor el modelo se ajusta mejor.

  • R-Squared: es otro indicador que informa de la bondad del ajuste del modelo a los datos. Ocila entre 0 y 1; valores próximos a 1 indican un buen ajuste del modelo lineal a los datos.

  • Adj R-Squared: similar a R², pero penaliza la introducción en el modelo de variables independientes poco relevantes a la hora de explicar la variable dependiente. Si R² = 0.9293 y R² ajustado = 0.9268, se puede concluir que el modelo lineal se ajusta de forma aceptable a nuestros datos

  • F-statistic: cuanto más alto mejor.

Merece la pena detenernos un poco más en los valores de los coeficientes. La tabla contiene:

  • Estimate. Como se ha señalado, son los coeficientes del modelo de regresión. Si están significativamente asociadas a la variable de resultado, están marcadas con estrellas.

  • Standard errors (SE): definen la precisión de los coeficientes anteriores. Para un coeficiente dado, el SE refleja cómo varía el coeficiente bajo muestreo repetido. Se puede utilizar para calcular los intervalos de confianza y la estadística t.

  • t-value: su valor se calcula dividiendo el coeficiente por la desviación típica. Es usado para verificar si el coeficiente es o no significativamente diferente de 0. Para ello, las hipótesis son:

    • Hipótesis Nula (H0): los coeficientes son iguales a 0 (no hay relación entre x e y).

    • Hipótesis alternativa (H1): los coeficientes no son iguales a 0 (si hay relación entre x e y).

Si no es significativamente diferente, el coeficiente no añade información al modelo y podríamos prescindir de él.

  • Pr(>|t|) es el nivel de significación. En nuestro ejemplo, ambos son menores que 0.05, luego considerando un nivel del significación alpha = 0.05, rechazamos la hipótesis nula en ambos contrastes, de manera que podemos suponer que ambas variables estan significativamente relacionadas.

4.4 Representación gráfica

Es habitual acompañar el modelo de regresión con un diagrama de puntos y dibujar sobre él la recta de regresión.

plot(mpg ~ wt, 
     data = mtcars, pch = 19, col = "darkblue")                      # Creamos el gráfico
abline(modelo_regresion_1, col = "red", lwd = 3)                                # Línea de regresión
text(paste("Correlación:", round(cor(mtcars$mpg, mtcars$wt), 2)), 
     x = 30,                                                                    # Posición etiqueta: sobre el valor 30 de x
     y = 5)                                                                     # Posición etiqueta: sobre el valor 5 de y

Incluso, sobre ese gráfico podemos incluir los coeficientes y el modelo de regresión

plot(mpg ~ wt, 
     data = mtcars, pch = 16)
abline(modelo_regresion_1, col="red", lwd=2)                                    # Añade una recta de regresión 
segments(mtcars$wt,                                                            # Añade las desviaciones individuales 
         fitted(modelo_regresion_1), 
         mtcars$wt, 
         mtcars$mpg, 
         col="red")

etiqueta_corr <- cor(mtcars$mpg,mtcars$wt)                                      # Calculamos la correlación entre las variables
text(paste("Correlación:", round(etiqueta_corr, 2)), 
     x = 5, 
     y = 30) 
coef <- round(coef(lm(mtcars$mpg ~ mtcars$wt)), 2)                              # Calculamos el modelo de regresion
text(paste("Y = ", coef[1], "+", coef[2], "x"),
     x = 5, 
     y = 27)

También es habitual dibujar una curva suavizada por medio del método LOWESS, para comparar sus resultados con los del método por mínimos cuadrados. Este método LOWESS es una regresión no paramética de carácter polinómico

plot(mpg ~ wt, data = mtcars, pch = 19, col = "darkblue") 
abline(modelo_regresion_1, col="red", lwd=2)                                    # Añade una recta de regresión 

lines(lowess(mtcars$wt, mtcars$mpg), col = "blue", lwd = 3)                     # Ajuste suavizado
legend("topright",                                                              
       legend = c("Lineal", "Suavizado"),
       lwd = 3, 
       lty = c(1, 1), 
       col = c("red", "blue"))

Desafío.

EJERCICIO 1. A partir del dataframe alumnos_UC.Rdata elabora un modelo de regresión que relacione las variables peso y altura.

  • Elabora sendos histogramas para evaluar la normalidad de la distribución de ambas variables.
media_edad <- mean(df$edad)                             # Cálculo de la media
desv_tipica_edad <- sd(df$edad)                         # Cálculo de la desviación típica
hist(df$edad, prob=TRUE,                               # Histograma
     xlab="edad", ylim=c(0, 0.30),
     main="Curva normal sobre el histograma")
curve(dnorm(x, mean=media_edad, sd=desv_tipica_edad),       # Curva normal
      col="darkblue", lwd=2, add=TRUE)

media_peso <- mean(df$peso)                             # Cálculo de la media
desv_tipica_peso <- sd(df$peso)                         # Cálculo de la desviación típica
hist(df$peso, prob=TRUE,                               # Histograma
     xlab="peso", ylim=c(0, 0.05),
     main="Curva normal sobre el histograma")
curve(dnorm(x, mean=media_peso, sd=desv_tipica_peso),       # Curva normal
      col="red", lwd=2, add=TRUE)
  • Somete ambas variables a un test de Shapiro para confirmar la normalidad de ambas variables.
shapiro.test(df$edad)          
shapiro.test(df$peso)
  • Visualiza en pantalla los valores ajustados (también llamados predicciones).
regresion_peso_edad <- lm(edad ~ peso,
                          data = df)

regresion_peso_edad$fitted.values

residuos_peso_edad <- regresion_peso_edad$fitted.values
df <- cbind(df, residuos_peso_edad)
  • A partir de la salida del resumen, obtén los valores de \(R^2\) y \(R^2 ajustado\)
summary(regresion_peso_edad)
  • Representa gráficamente la relación entre ambas variables con un diagrama de puntos sobre el que deberás dibujar la recta de regresión, los coeficientes y el modelo de regresión
plot(edad ~ peso,
     data = df, pch = 16)
abline(regresion_peso_edad, col="red", lwd=2)                                    # Añade una recta de regresión

etiqueta_corr <- cor(df$peso,df$edad)                                      # Calculamos la correlación entre las variables
text(paste("Correlación:", round(etiqueta_corr, 2)),
     x = 60,
     y = 24)
coef <- round(coef(lm(df$edad ~ df$peso)), 2)                              # Calculamos el modelo de regresion
text(paste("Y = ", coef[1], "+", coef[2], "x"),
     x = 70,
     y = 20)

EJERCICIO 2. A partir del dataframe zonas_verdes.Rdata elabora un modelo de regresión que relacione las variables superficie (variable dependiente) y parques (variable independiente).

  • Elabora un gráfico q-q para evaluar la normalidad de la distribución de ambas variables.
library("car")
qqPlot(zonas_verdes$superficie)              
qqPlot(zonas_verdes$parques)
  • Somete ambas variables a un test de Shapiro para confirmar la normalidad de ambas variables.
shapiro.test(zonas_verdes$superficie)          
shapiro.test(zonas_verdes$parques)
  • Visualiza en pantalla los resíduos del modelo.
regresion_superficie_parques <- lm(superficie ~ parques,
                          data = zonas_verdes)

observados_vs_prediccion <- cbind(regresion_superficie_parques$residuals, regresion_superficie_parques$fitted, zonas_verdes$parques)
  • Representa gráficamente la relación entre ambas variables con un diagrama de puntos sobre el que deberás dibujar la recta de regresión, los coeficientes y el modelo de regresión
plot(superficie ~ parques,
     data = zonas_verdes, pch = 10)
abline(regresion_superficie_parques, col="red", lwd=2)

4.5 Estadísticos y gráficos de diagnóstico del modelo

Que observemos cómo la recta dibujada anteriormente parece ajustarse a los datos iniciales no significa que el modelo sea correcto:

  • Si sólo pretendemos hallar la relación entre dos variables, con calcular el modelo ya sería suficiente.

  • Si el fin del modelo es predecir, es necesario comprobar que se verifican unas reglas ya establecidas que aseguran la bondad (calidad) del modelo.

Esta comprobación se puede realizar de dos maneras:

  • Calculando una serie de estadísticos.

  • Representando gráficamente esos estadísticos mediante una serie de gráficos de diagnóstico. Estos gráficos se aplican a los resíduos del modelo de regresión, e identifican posibles problemas en el modelo que podrían en riesgo los resultados obtenidos. Esto se hace con la función plot() aplicada al objeto resultante de aplicar la función lm().

Podemos comenzar representando estos gráficos.

plot(modelo_regresion_1)

Estos gráficos también se pueden obtener en un único panel:

par(mfrow=c(2,2))                       
plot(modelo_regresion_1, which = 1)
plot(modelo_regresion_1, which = 2)
plot(modelo_regresion_1, which = 3)
plot(modelo_regresion_1, which = 4)
par(mfrow=c(1,1))                       # devuelve la pantalla a su estado original

El significado de estos gráficos es el siguiente:

  • Gráfico 1: RESIDUALS vs. LEVERAGE. Sirve para identificar observaciones influyentes, que son aquellas que aparecen fuera de la línea discontínua, que a su vez reproduce la distancia de Cook. Detecta puntos con una influencia importante en el cálculo de las estimaciones de los parámetros. En caso de aparecer estas observaciones, debe verificarse su validez, planteándose su eliminación o la búsqueda de observaciones similares para incorporar en el modelo. En el ejemplo, el modelo Chrysler Imperial se encuentra más allá de esa línea, lo que lo convierte en una observación influyente.

Podemos mejorar la identificación de casos inusuales (gráfico 1) recurriendo a la representación de un boxplot.

par(mfrow=c(1, 2))                                      # Dividimos el área de los gráficos en 2 columnas
boxplot(mtcars$mpg, main="Variable X")                  # Gráfico de caja y bigotes para la variable X
boxplot(mtcars$wt, main="Variable Y")

par(mfrow=c(1,1))                                      # Devolvemos la pantalla a su estado original

Valores atípicos - La independencia de los errores

residuos <- residuals(modelo_regresion_1)                                       # Residuos
plot(residuos,
     rstandard(modelo_regresion_1),
     xlab="Variable X",
     ylab="Variable Y")

Puede observarse que, de acuerdo a la definición de casos atípicos, la variable wtcontaría con dos observaciones de este tipo. Para confirmar que ambos casos si pueden ser considerados atípicos, los cuantificaremos con valores elevados de parámetro Studentised residual obtenida con la función rstudent().

outliers <- rstudent(modelo_regresion_1)
outliers

También podemos obtener el valor del Leverage, que señala a observaciones muy diferentes de las restantes, aunque no corresponde necesariamente a un residuo muy alto.

leverage <- hatvalues(modelo_regresion_1)
leverage

Y finalmente, la distancia de Cook.

distancia_cook <- cooks.distance(modelo_regresion_1)
distancia_cook

Valores superiores a 1 en la distancia de Cook’son considerados muy grandes

plot(cooks.distance(modelo_regresion_1), 
     pch = 16, 
     col = "blue")
  • Gráfico 2: SCALE vs LOCATION

Este gráfico se usa para verificar la asumción de igualdad de varianzas (también conocida como “homoscedasticidad”) entre los resíduos del modelo de regresión. Si la línea roja que aparece en el gráfico es aproximadamente horizontal, se cumple esta condición; idealmente, los residuos deben estar aleatoriamente distribuidos a lo largo del gráfico, sin formar ningún tipo de patrón.

Figura 3: Gráfico de los resíduos de un modelo de regresión mostrando heteroscedasticidad

En el ejemplo citado se ve cómo la línea no es totalmente horizontal, si bien no se desvía demasiado de la horizontalidad. Entonces se puede considerar que el modelo no viola la asumción de igualdad de varianzas.

Es posible desarrollar un poco más esta condición:

plot(fitted.values(modelo_regresion_1),
     rstandard(modelo_regresion_1), 
     xlab="Valores ajustados", 
     ylab="Residuos estandarizados",
     ylim = c(-3,3))                                        # gráfico 2D de los valores ajustados vs. los residuos estandarizados 
abline(h=0, col = "red")

Además, se puede aplicar la prueba de Breush-Pagan. Si el valor-p es menor que el nivel de significancia usual de 5%, no se cumple la homocedasticidad.

library(lmtest)                                           
bptest(modelo_regresion_1)

Una aproximación aún más correcta es someter los residuos al test de varianza no constante (Cook y Weisberg 1983). De acuerdo con esta prueba:

  • H0 hipótesis nula: los residuos violan la homogeneidad de la varianza.

  • H1 Hipotesis alternativa: los residuos no violan la homogeneidad de la varianza.

Probabilidades superiores a 0,05 rechazan la hipótesis nula: son indicativo de que los residuos no violan la homogeneidad de la varianza

if(!require("car")) install.packages("car")
library(car)
residualPlots(modelo_regresion_1) 
ncvTest(modelo_regresion_1)
  • Gráfico 3: gráfico Q-Q. Este gráfico se usa para determinar si los resíduos del modelo de regresión poseen una distribución normal. Si los puntos, que representan los casos, se sitúan sobre la recta, se asume que los resíduos están distribuidos normalmente.Los puntos alejados de la recta muestran una etiqueta o su número de caso (filas del dataframe). En nuestro modelo la mayoría de los puntos se superponen a la diagonal, si bien algunas observaciones (Fiat 128, Torota Corolla, Chrysler Imperial) se alejan de la línea en la terminación. Esto no impide que consideremos que los residuos tienen una distribución normal.

Otras posibilidades de evaluación de la normalidad de los residuos son las siguientes:

hist(residuos,                          
     xlab = "Valor del resíduo",        
     main = "Histograma de resíduos",                       
     breaks = 10)                      

qqnorm(residuos)
qqline(residuos)

Podrían también transformarse en puntuaciones z. Es deseable que los residuos estandarizados estén lo más cerca posible a la línea punteada que aparece en el gráfico.

z_residuos <- rstandard(modelo_regresion_1)         # Calcula los resíduos estandarizados del modelo ajustado (completo) 

par(mfrow=c(1,3))                          # Divide la ventana en una fila y tres columnas 
hist(z_residuos)                           # Histograma de los residuos estandarizados 
boxplot(z_residuos, ylim = c(-3,3))        # Diagrama de cajas de los residuos estandarizados 
qqnorm(z_residuos)                         # gráfico de cuantiles de los residuos estandarizados
qqline(z_residuos)                         

par(mfrow=c(1,1))                       # devuelve la pantalla a su estado original

Y por último, ser sometidos a un contraste de hipótesis.

shapiro.test(modelo_regresion_1$residuals)

De acuerdo con todo lo expuesto, los resíduos presentan una distribución que no se aleja de la distribución normal.

  • Gráfico 4: RESIDUALS vs. FITTED (residuos vs valores ajustados) Este gráfico sirve para determinar si los resíduos muestran algún tipo de patrón no lineal. Si la línea roja atraviesa horizontalmente el gráfico se asume que los resíduos siguen un patrón lineal y no muestran ningún patrón sistemático (por ejemplo, un arco) que señalaría problemas en el modelo (por ejemplo, una nube de puntos). En el modelo obtenido, la línea roja (filtro lowess) se desvía bastante de una disposición horizontal, por lo que habría que analizar si un modelo de regresión lineal por minimos cuadrados es el método más conveniente para estos datos.

Linearidad de la relación

valores_ajustados_1 <- fitted.values(modelo_regresion_1)
plot(valores_ajustados_1, 
     mtcars$wt,
     xlab = "Fitted Values",
     ylab = "Observed Values")

Si se quisiera realizar un análisis en profundidad de todos estos aspectos, se podrían integrar en un único dataframe.

datos_problema_regresion <- data.frame(mtcars$mpg,
                                        mtcars$wt,
                                        modelo_regresion_1$fitted.values,
                                        modelo_regresion_1$residuals,
                                        outliers,
                                        leverage,
                                        distancia_cook)

4.6 Obtención de predicciones

La función predict() se usa para predecir el valor de la variable dependiente para nuevas observaciones de la variable independiente. Para ello, se creará un nuevo dataframe que contendrá 10 valores nuevos (que oscilan entre 1.5 y 5,5) de la variable independiente wt. Queremos saber cuánto podrán circular estos nuevos vehículos utilizando el modelo de regresión calculado previamente.

nuevos_datos <- data.frame(wt = sample(1:5.5, size=10, replace = TRUE))
nuevos_datos

Estos nuevos datos van a incorporarse a la función predict(), que usa la información proporcionada por el modelo de regresión para predecir el valor teórico de \(y\) (mpg) en función de los nuevos valores de \(x\) (wt). Al mismo tiempo, se creará un nuevo objeto que incluirá los nuevos valores de \(x\) con los valores teóricos de \(y\).

wt_predicciones <- predict(modelo_regresion_1,                                               # Nuevo dataframe con valores nuevos de y 
                     nuevos_datos)   
nuevos_datos <- cbind(nuevos_datos, wt_predicciones)

5 INFORMACIÓN ADICIONAL

La evaluación de un modelo de regresión también puede realizarse sobre un nuevo conjunto de datos que no ha sido incluido en el desarrollo del modelo. Para ilustrar este procedimiento incorporamos el dataset Boston del paquete MASS. Este paquete incluye la variable age, que corresponde a la proporción de viviendas en propiedad construidas antes de 1940.

library("MASS")
data("Boston")                # Cargamos el fichero directamente del conjunto de ficheros disponibles en R
View(Boston)

Para determinar qué variables están mejor correlacionadas con age, calculamos una matriz de correlaciones entre todas las variables.

round(cor(Boston), 2)

Los resultados muestran que age está correlacionada positivamente con nox (contaminante relacionado con el tráfico de vehículos) y negativamente con dis (distancia promedio a los grandes distritos comerciales). En este último caso, podemos suponer que la mayoría de las viviendas más antiguas se sitúan en el centro de las ciudades, alejadas de los actuales centros de negocios, que suelen situarse en la periferia.

En consecuencia, crearemos un nuevo dataframe sólo con ambas variables.

boston <- subset(Boston, select = c(age, dis))
rm(Boston)

El procedimiento habitual el siguiente:

?floor
smp_size <- floor(0.75 * nrow(boston))    # Devuelve el número entero más pequeño que esa fórmula: conjunto de entrenamiento
set.seed(123)                                                     # Fija la semilla para poder reproducir más tarde el procedimiento
train_ind <- sample(seq_len(nrow(boston)), 
                    size = smp_size)
entrenamiento.1 <- boston[train_ind,     ]
verificacion.1 <- boston[-train_ind,     ]

Otra manera de hacerlo es mediante la creación de una variable de apoyo para crear los dos subconjuntos (70% vs 30%)

variable_apoyo <- sample(c(rep(0, 0.7 * nrow(boston)),            
                           rep(1, 0.3 * nrow(boston))))
variable_apoyo                                                                  # Verificar el contenido de la variable de apoyo
table(variable_apoyo)                                                           # Verificar si las frecuencias son correctas   
entrenamiento.2 <- boston[variable_apoyo == 0, ]                                # Subconjunto de entrenamiento
verificacion.2 <- boston[variable_apoyo == 1, ]                             # Subconjunto de verificación
modelo.regresion.2 <- lm(age ~ dis, 
                         data = entrenamiento.1)

Resumen del modelo de regresión

summary(modelo.regresion.2)

Predicciones

library(caret)
prediccion.dis <- predict(modelo.regresion.2,
                         verificacion.1)   
verificacion.1 <- cbind(verificacion.1, 
                        prediccion.dis)

Evaluación del modelo en los datos de verificación. La raíz del error cuadrático medio (RSME): corresponde a la diferencia media entre los valores predichos y los valores observados. Cuanto más bajo se este error, mejor es el modelo.

RMSE(verificacion.1$dis, verificacion.1$prediccion.dis)         
R2(verificacion.1$dis, verificacion.1$prediccion.dis)     # Correlación entre los valores observados y los predichos (cuanto más altos, mejor).

ACTIVIDADES DE EVALUACIÓN CONTINUA

El fichero word 03_Actividad_evaluacion_continua_3_2023 contiene una serie de actividades encaminadas a poner en práctica los conocimientos adquiridos. Cada alumnos deberá elaborar un script que contenga los procedimientos necesarios para responder a las cuestiones planteadas, y enviarlo posteriormente al profesor a través de correo electrónico.

Los resultados de esta actividad pueden consultarse en 03_Actividad_evaluacion_continua_3_2023_resuelta