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:
\(a\): es la pendiente de la recta. Indica cuánto aumenta \(y\) por cada aumento de una unidad de \(x\).
\(b\): es el origen de la recta. Equivale al valor de \(y\) cuando \(x\) es 0.
\(e\): es el error aleatorio (también conocidos como residuos). Son las diferencias entre el valor teórico proporcionado por el modelo estadístico y el valor real (equivalente a aquella parte de la relación que no puede ser explicada por el modelo), y se atribuye a errores de medición o desviaciones aleatorias de las condiciones de regresión.
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.
Supongamos las siguientes variables:
x_ejemplo1 <- c(0.8,1,1.2,1.3)
y_ejemplo1 <-c(1,2,3,5)
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
## [1] 7.288136
# Parámetro B
B_ejemplo1 <- media_y_ejemplo1 - A_ejemplo1 * media_x_ejemplo1
B_ejemplo1
## [1] -5.084746
Por lo tanto, nuestro modelo de regresión tomaría la forma de \(y = 7.2x - 5.084746\).
rm(x_ejemplo1, y_ejemplo1)
La función básica para obtener un modelo de regresión lineal es
lm()
. Su sintaxis es la siguiente:
lm(formula, data, …)
donde:
y ~ x
en la que se
especifica cuál es la variable dependiente (\(y\)) y cuál es la variable independiente
(\(x\)).data
especifica de qué fichero se
extraen dichas variables.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)
##
## Call:
## lm(formula = y_ejemplo1 ~ x_ejemplo1)
##
## Coefficients:
## (Intercept) x_ejemplo1
## -5.085 7.288
En la salida del procedimiento, los coefficents
equivalen a los parámetros a
y b
.
b = -5.085. Es el origen. Este valor equivale al valor de la variable y cuando el valor de la variable x es 0. En muchas ocasiones, este parámetro no es demasiado útil, dado que un valor \(x = 0\) es absolutamente teórico.
a = 7.288. Es la pendiente. Esto significa que un aumento de una unidad en \(x\) está asociado con un aumento de 7,288 unidades en \(y\).
Aunque parezca sencillo, la elaboración de un modelo de regresión requiere superar varios pasos:
Evaluar si las variables \(x\) e \(y\) cumplen ciertos supuestos.
Obtener el modelo.
Analizar los diferentes parámetros que resumen el modelo de regresión.
Revisar ciertos estadísticos y gráficos de diagnóstico.
Dibujar el modelo de regresión.
Utilizar el modelo para hacer predicciones.
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:
m
: consumo del vehículo (recorrido en millas con un
galón de combustible -aprox 3.5 litros).
wt
: peso del vehículo (1000 lbs)
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)
## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...
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
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)
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):
## [1] 20 18
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-Wilk normality test
##
## data: mtcars$mpg
## W = 0.94756, p-value = 0.1229
shapiro.test(mtcars$wt)
##
## Shapiro-Wilk normality test
##
## data: mtcars$wt
## W = 0.94326, p-value = 0.09265
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, desv_tipica_wt, media_mpg, media_wt)
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)
##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Coefficients:
## (Intercept) wt
## 37.285 -5.344
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)
## $names
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
##
## $class
## [1] "lm"
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
## (Intercept) wt
## 37.285126 -5.344472
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
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## 23.282611 21.919770 24.885952 20.102650
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 18.900144 18.793255 18.205363 20.236262
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 20.450041 18.900144 18.900144 15.533127
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 17.350247 17.083024 9.226650 8.296712
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 8.718926 25.527289 28.653805 27.478021
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## 24.111004 18.472586 18.926866 16.762355
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 16.735633 26.943574 25.847957 29.198941
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## 20.343151 22.480940 18.205363 22.427495
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
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## -2.2826106 -0.9197704 -2.0859521 1.2973499
## Hornet Sportabout Valiant Duster 360 Merc 240D
## -0.2001440 -0.6932545 -3.9053627 4.1637381
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 2.3499593 0.2998560 -1.1001440 0.8668731
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## -0.0502472 -1.8830236 1.1733496 2.1032876
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 5.9810744 6.8727113 1.7461954 6.4219792
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## -2.6110037 -2.9725862 -3.7268663 -3.4623553
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 2.4643670 0.3564263 0.1520430 1.2010593
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## -4.5431513 -2.7809399 -3.2053627 -1.0274952
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)
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## 23.282611 21.919770 24.885952 20.102650
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 18.900144 18.793255 18.205363 20.236262
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 20.450041 18.900144 18.900144 15.533127
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 17.350247 17.083024 9.226650 8.296712
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 8.718926 25.527289 28.653805 27.478021
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## 24.111004 18.472586 18.926866 16.762355
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 16.735633 26.943574 25.847957 29.198941
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## 20.343151 22.480940 18.205363 22.427495
residuals(modelo_regresion_1)
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## -2.2826106 -0.9197704 -2.0859521 1.2973499
## Hornet Sportabout Valiant Duster 360 Merc 240D
## -0.2001440 -0.6932545 -3.9053627 4.1637381
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 2.3499593 0.2998560 -1.1001440 0.8668731
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## -0.0502472 -1.8830236 1.1733496 2.1032876
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 5.9810744 6.8727113 1.7461954 6.4219792
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## -2.6110037 -2.9725862 -3.7268663 -3.4623553
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 2.4643670 0.3564263 0.1520430 1.2010593
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## -4.5431513 -2.7809399 -3.2053627 -1.0274952
El resumen del modelo de regresión se obtiene aplicando la función
summary
.
summary(modelo_regresion_1)
##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10
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.
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"))
📝 ACTIVIDAD DE EVALUACIÓN CONTINUA: :
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.
Somete ambas variables a un test de Shapiro para confirmar la normalidad de ambas variables.
Visualiza en pantalla los valores ajustados (también llamados predicciones).
A partir de la salida del resumen, obtén los valores de \(R^2\) y \(R^2 ajustado\)
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
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.
Somete ambas variables a un test de Shapiro para confirmar la normalidad de ambas variables.
Visualiza en pantalla los resíduos del modelo.
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
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:
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 wt
contarí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
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## -0.76076977 -0.30274050 -0.69972329 0.42681176
## Hornet Sportabout Valiant Duster 360 Merc 240D
## -0.06570059 -0.22779587 -1.32167204 1.41169344
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 0.77876838 0.09844187 -0.36192838 0.28837628
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## -0.01655496 -0.62515414 0.41709709 0.76444774
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 2.32816206 2.53780106 0.60385911 2.38384376
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## -0.87622897 -0.99290217 -1.25610443 -1.16991451
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 0.82326228 0.12041613 0.05090396 0.41668241
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## -1.55056010 -0.92873434 -1.07426703 -0.33877056
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
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## 0.04326896 0.03519677 0.05837573 0.03125017
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 0.03292182 0.03323551 0.03544265 0.03127502
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 0.03140238 0.03292182 0.03292182 0.05575179
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 0.04010861 0.04192052 0.17047665 0.19533191
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 0.18379417 0.06611662 0.11774978 0.09562654
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## 0.05031684 0.03433832 0.03284761 0.04431718
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 0.04452785 0.08664873 0.07035096 0.12911356
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## 0.03132522 0.03798993 0.03544265 0.03769190
Y finalmente, la distancia de Cook.
distancia_cook <- cooks.distance(modelo_regresion_1)
distancia_cook
## Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive
## 1.327407e-02 1.723963e-03 1.543937e-02 3.020558e-03
## Hornet Sportabout Valiant Duster 360 Merc 240D
## 7.599578e-05 9.210650e-04 3.131386e-02 3.113918e-02
## Merc 230 Merc 280 Merc 280C Merc 450SE
## 9.961857e-03 1.705808e-04 2.296167e-03 2.532453e-03
## Merc 450SL Merc 450SLC Cadillac Fleetwood Lincoln Continental
## 5.923263e-06 8.727280e-03 1.838257e-02 7.192515e-02
## Chrysler Imperial Fiat 128 Honda Civic Toyota Corolla
## 5.319056e-01 1.929858e-01 2.486028e-02 2.598750e-01
## Toyota Corona Dodge Challenger AMC Javelin Camaro Z28
## 2.049815e-02 1.753646e-02 2.628727e-02 3.134958e-02
## Pontiac Firebird Fiat X1-9 Porsche 914-2 Lotus Europa
## 1.596432e-02 7.111643e-04 1.014165e-04 1.323493e-02
## Ford Pantera L Ferrari Dino Maserati Bora Volvo 142E
## 3.713611e-02 1.710946e-02 2.109445e-02 2.315918e-03
Valores superiores a 1 en la distancia de Cook’son considerados muy grandes
plot(cooks.distance(modelo_regresion_1),
pch = 16,
col = "blue")
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.
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)
##
## studentized Breusch-Pagan test
##
## data: modelo_regresion_1
## BP = 0.040438, df = 1, p-value = 0.8406
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)
## Test stat Pr(>|Test stat|)
## wt 3.258 0.002860 **
## Tukey test 3.258 0.001122 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ncvTest(modelo_regresion_1)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.03794177, Df = 1, p = 0.84556
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)
##
## Shapiro-Wilk normality test
##
## data: modelo_regresion_1$residuals
## W = 0.94508, p-value = 0.1044
De acuerdo con todo lo expuesto, los resíduos presentan una distribución que no se aleja de la distribución normal.
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)
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
## wt
## 1 4
## 2 2
## 3 5
## 4 4
## 5 1
## 6 5
## 7 4
## 8 2
## 9 2
## 10 4
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)
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)
## crim zn indus chas nox rm age dis rad tax ptratio
## crim 1.00 -0.20 0.41 -0.06 0.42 -0.22 0.35 -0.38 0.63 0.58 0.29
## zn -0.20 1.00 -0.53 -0.04 -0.52 0.31 -0.57 0.66 -0.31 -0.31 -0.39
## indus 0.41 -0.53 1.00 0.06 0.76 -0.39 0.64 -0.71 0.60 0.72 0.38
## chas -0.06 -0.04 0.06 1.00 0.09 0.09 0.09 -0.10 -0.01 -0.04 -0.12
## nox 0.42 -0.52 0.76 0.09 1.00 -0.30 0.73 -0.77 0.61 0.67 0.19
## rm -0.22 0.31 -0.39 0.09 -0.30 1.00 -0.24 0.21 -0.21 -0.29 -0.36
## age 0.35 -0.57 0.64 0.09 0.73 -0.24 1.00 -0.75 0.46 0.51 0.26
## dis -0.38 0.66 -0.71 -0.10 -0.77 0.21 -0.75 1.00 -0.49 -0.53 -0.23
## rad 0.63 -0.31 0.60 -0.01 0.61 -0.21 0.46 -0.49 1.00 0.91 0.46
## tax 0.58 -0.31 0.72 -0.04 0.67 -0.29 0.51 -0.53 0.91 1.00 0.46
## ptratio 0.29 -0.39 0.38 -0.12 0.19 -0.36 0.26 -0.23 0.46 0.46 1.00
## black -0.39 0.18 -0.36 0.05 -0.38 0.13 -0.27 0.29 -0.44 -0.44 -0.18
## lstat 0.46 -0.41 0.60 -0.05 0.59 -0.61 0.60 -0.50 0.49 0.54 0.37
## medv -0.39 0.36 -0.48 0.18 -0.43 0.70 -0.38 0.25 -0.38 -0.47 -0.51
## black lstat medv
## crim -0.39 0.46 -0.39
## zn 0.18 -0.41 0.36
## indus -0.36 0.60 -0.48
## chas 0.05 -0.05 0.18
## nox -0.38 0.59 -0.43
## rm 0.13 -0.61 0.70
## age -0.27 0.60 -0.38
## dis 0.29 -0.50 0.25
## rad -0.44 0.49 -0.38
## tax -0.44 0.54 -0.47
## ptratio -0.18 0.37 -0.51
## black 1.00 -0.37 0.33
## lstat -0.37 1.00 -0.74
## medv 0.33 -0.74 1.00
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
## [1] 0 0 1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1
## [38] 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 1 0 0
## [75] 0 0 1 0 0 1 0 0 0 1 1 0 1 0 0 0 1 0 1 1 0 1 1 0 0 0 0 0 1 0 0 0 1 1 1 0 0
## [112] 1 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 1 1 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 1 1 0 0
## [149] 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 1
## [186] 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1
## [223] 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 0 1 0 0 0
## [260] 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 1 0 1 0 1 0 1 1 1 1 0 1 1 0 0 0 0 0 1 0 0 0
## [297] 1 0 0 1 0 0 0 1 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 1 0 0 0
## [334] 0 0 1 0 1 1 1 1 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 1 1
## [371] 1 0 1 0 0 0 1 1 0 0 1 1 0 0 0 0 1 0 0 0 1 1 0 1 0 0 1 1 0 1 0 0 0 0 0 0 0
## [408] 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 1 1 1 0 1 0
## [445] 0 0 1 1 0 0 1 0 0 0 0 0 1 1 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## [482] 0 0 1 0 1 0 0 0 0 1 0 0 1 1 1 1 0 0 0 0 1 0 1 0
table(variable_apoyo) # Verificar si las frecuencias son correctas
## variable_apoyo
## 0 1
## 354 151
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)
##
## Call:
## lm(formula = age ~ dis, data = entrenamiento.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.960 -9.377 3.993 9.685 53.858
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 107.2272 1.9255 55.69 <2e-16 ***
## dis -10.1651 0.4362 -23.30 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.27 on 377 degrees of freedom
## Multiple R-squared: 0.5903, Adjusted R-squared: 0.5892
## F-statistic: 543.1 on 1 and 377 DF, p-value: < 2.2e-16
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)
## [1] 70.27605
R2(verificacion.1$dis, verificacion.1$prediccion.dis) # Correlación entre los valores observados y los predichos (cuanto más altos, mejor).
## [1] 1
📝 ACTIVIDAD DE EVALUACIÓN CONTINUA: :
EJERCICIO 3. A partir de la hoja Ejemplo2
presente en la
hoja de cálculo https://personales.unican.es/rasillad/docencia/G14/evaluacion_continua/ejercicios_regresion_2022.xlsx,
realiza los siguientes cálculos:
library(readxl)
download.file("https://personales.unican.es/rasillad/docencia/G14/TEMA_3/ejercicios_regresion_2023.xlsx",
destfile = "regresion.xlsx",
mode="wb")
ejercicio2 <- read_excel("regresion.xlsx", sheet = "Ejemplo2")
ejercicio2
## # A tibble: 25 × 3
## peso edad grasas
## <dbl> <dbl> <dbl>
## 1 84 46 354
## 2 73 20 190
## 3 65 52 405
## 4 70 30 263
## 5 76 57 451
## 6 69 25 302
## 7 63 28 288
## 8 72 36 385
## 9 79 57 402
## 10 75 44 365
## # ℹ 15 more rows
Representa gráficamente la relación entre cada una de las variables (peso, edad y grasas) con las restantes, para determinar si existe una relación lineal entre ellas. La relación lineal podemos observarla mediante un gráfico
Comprueba si esas variables cumplen en supuesto de normalidad.
Teniendo en cuenta que el contenido en grasas puede considerarse dependiente de la edad y el peso, identifica cuál está mejor correlacionada con las grasas y desarrolla un modelo de regresión lineal.
EJERCICIO 4. A partir de la hoja difusión presente en la hoja de cálculo https://personales.unican.es/rasillad/docencia/G14/TEMA_3/ejercicios_regresion_2023.xlsx, realiza los siguientes cálculos:
Crea un data frame que almacene los datos del ejemplo propuesto.
Representa gráficamente la relación entre la variable ventas y las restantes variables, para determinar si existe una relación lineal entre ellas.
Comprueba si todas las variables cumplen en supuesto de normalidad.
Calcula una matriz de correlación entre esas variables y represéntalas gráficamente mediante un correlograma.
Si la variable ventas puede considerarse la variable dependiente, identifica con qué variable está mejor correlacionada y desarrolla un modelo de regresión lineal.
Interpreta los parámetros de esa recta ¿son significativos estos parámetros? ¿Qué puede decirse del ajuste del modelo a los datos?
Si para el mes que viene quisiéramos elevar nuestras ventas a 25 millones de euros, ¿cuál sería el valor de la variable independiente que nos permitiría alcanzar este valor?
EJERCICIO 5. A partir de la hoja clima_madrid
presente
en la hoja de cálculo https://personales.unican.es/rasillad/docencia/G14/TEMA_3/ejercicios_regresion_2023.xlsx,
realiza los siguientes cálculos:
Crea un data frame que almacene los datos del ejemplo propuesto.
Representa gráficamente la relación entre las variables alt, pp, tene y tjul, para determinar si existe una relación lineal entre ellas.
Comprueba si esas variables cumplen en supuesto de normalidad. Existen diferentes posibilidades. Una de ellas es el gráfico Q-Q.
Otra posibilidad es el test de Shapiro-Wilks
Calcula una matriz de correlación entre todas las variables y elabora al menos dos gráficos que te permitan observar de manera global el grado de relación entre todas esas variables.
Teniendo en cuenta que la altitud es una variable geográfica que genera alteraciones importantes en la distribución espacial de las variables climáticas, identifica cuál está mejor correlacionada con la altitud y desarrolla un modelo de regresión lineal.
Interpreta los parámetros de esa recta ¿son significativos estos parámetros? ¿Qué puede decirse del ajuste del modelo a los datos?
Queremos trasladarnos a una vivienda que estará situada en el Puerto de Somosierra a 1480 metros de altitud. ¿Qué precipitación registraremos?
EJERCICIO 6. El dataset “Air quality” contiene los valores diarios de las siguientes variables ambientales entre el 1 de mayo y el 30 de septiembre de 1973. Se puede acceder al mismo desde R con los siguientes comandos:
Las variables de ese conjunto de datos son las siguientes:
• Month: mes.
• Day: día.
• Ozone: media diaria del ozono (en partes por millón) entre las 13 y las 15 horas
• Solar.R: radiación solar (en Langleys) entre las 08 y las 12 horas.
• Wind: velocidad media del viento (en millas por hora) entre las 07 y las 10.
• Temp: temperatura máxima diaria en ºF.
Transforma los valores de la temperatura máxima y viento en unidades del sistema métrico decimal. Para ello:
Aplica a la temperatura del aire el siguiente factor de conversión: temperatura en ºF multiplicada por 0,5556 – 32.
Aplica al viento el siguiente factor de conversión: viento en millas * 0,44704.
Convierte en factores las variables día y mes.
Realiza un análisis exploratorio de las variables ozono, radiación.solar, viento y temperatura.máxima en el que analiza cómo varían esos valores a lo largo de los meses analizados. Utiliza las herramientas que consideres más oportunas. Comenta brevemente los resultados.
Analiza si las variables ozono, radiación.solar, viento y temperatura máxima tienen una distribución normal por los métodos que consideres oportunos o si presentan algún tipo de sesgo. Si existiera alguna variable que no sigue una distribución normal, transfórmala elevando sus valores originales al cuadrado y repite el análisis para comprobar si la distribución ha cambiado.
Identifica que parejas de variables poseen un valor del coeficiente de correlación.
Teniendo en cuenta que consideramos la variable ozono como variable dependiente, elabora un diagrama de puntos que muestre la relación entre el ozono y la variable independiente mejor correlacionada.
Elabora un modelo de regresión entre el ozono como variable dependiente y la variable independiente mejor correlacionada.