Un modelo de regresión lineal múltiple para la predicción de ventas - la violación de la normalidad
Introducción al problema
Vamos a crear un modelo de regresión lineal que prediga las ventas en función del gasto en diferentes conceptos. La idea de este modelo es mostrar qué sucede cuando ignoramos uno de los requisitos fundamentales de la regresión lineal: la normalidad de los residuos.
Veremos que, pese a todo, somos capaces de construir un modelo predictivo decente que nos puede ayudar a la comprensión del problema.
Los datos pueden descargarse de aquí.
Asumiendo que ya tenemos los datos, vamos a cargarlos y a empezar por echar un vistazo a posibles outliers, valores sin informar, …
#cargamos los datos
advertising <- read.csv("advertising.csv")
head(advertising,6)## TV Radio Newspaper Sales
## 1 230.1 37.8 69.2 22.1
## 2 44.5 39.3 45.1 10.4
## 3 17.2 45.9 69.3 12.0
## 4 151.5 41.3 58.5 16.5
## 5 180.8 10.8 58.4 17.9
## 6 8.7 48.9 75.0 7.2Análisis exploratorio preliminar de los datos
Vamos a ver si hay algún N/A o outlier.
summary(advertising)## TV Radio Newspaper Sales
## Min. : 0.70 Min. : 0.000 Min. : 0.30 Min. : 1.60
## 1st Qu.: 74.38 1st Qu.: 9.975 1st Qu.: 12.75 1st Qu.:11.00
## Median :149.75 Median :22.900 Median : 25.75 Median :16.00
## Mean :147.04 Mean :23.264 Mean : 30.55 Mean :15.13
## 3rd Qu.:218.82 3rd Qu.:36.525 3rd Qu.: 45.10 3rd Qu.:19.05
## Max. :296.40 Max. :49.600 Max. :114.00 Max. :27.00En principio los datos están bien. No tenemos NA y todas nuestras variables son numéricas. Nada mejor que echarle un vistazo a los boxplot para buscar posibles outliers:
par(mfrow=c(1,3))
boxplot(advertising$TV,main="gasto en TV")
boxplot(advertising$Radio, main="gasto en Radio")
boxplot(advertising$Newspaper, main="gasto en Periodicos")
De la tercera figura, observamos que hay dos outliers en Newspaper. Son,
como vemos a continuación, 114 y 100.9. No tendría sentido hacer nada
con ellos, dado que es posible que sean valores legítimos y no tenemos
más información al respecto.
head(sort(advertising$Newspaper,decreasing=TRUE))## [1] 114.0 100.9 89.4 84.8 79.2 75.6Correlación lineal de cada una de las variables con las ventas
Antes de lanzarnos a crear un modelo matemático para el problema, es necesario comprobar que existe una relación lineal sólida entre cada una de las variables predictoras y la variable de salida, Sales.
library(corrplot)
M<-cor(advertising)
corrplot(M, method="circle")
Está claro que la mayor dependencia es con TV. Iremos quitando una a una las demás hasta conseguir que todos los coeficientes de regresión sean significativos.
Aparentemente, no tenemos problemas de multicolinealidad, ya que la correlación entre las distintas variables regresoras es pequeña.
Diagramas de dispersión para cada variable junto con Sales
A la vista del gráfico de correlaciones de arriba, tendremos una fuerte correlación lineal de las ventas con TV, más debil con Radio y practicamente nula con Newspaper. Veamos esto de forma gráfica:
par(mfrow=c(2,2))
plot(advertising$TV,advertising$Sales)
plot(advertising$Radio,advertising$Sales)
plot(advertising$Newspaper,advertising$Sales)
Puede observarse que la varianza no es constante en el primer gráfico.
Así como una correlación mucho más débil en los otros dos de Sales con
respecto de Radio y Newspaper.
Separando datos de test y de validación
En cada uno de los modelos de regresión, en general machine learning - deberemos siempre separar dos particiones. Una para entrenar el modelo y otra para validarlo una vez calculado.
Lo que hacemos aquí debajo es separar los datos en dos partes, train y test. La primera para entrenar y la segunda para validar.
#make this example reproducible
library(caTools)
set.seed(123)
#usamos el 90% de los datos para crear el modelo y el 20% para su validación.
sample <- sample.split(advertising$Sales, SplitRatio = 0.8)
train <- subset(advertising, sample == TRUE)
test <- subset(advertising, sample == FALSE)
nrow(train) ## [1] 160nrow(test)## [1] 40En este caso, hemos separado los datos en una parte de entrenamiento, que es aproximadamente el 80% de los datos iniciales y una parte de validación, que son los datos restantes. Los modelos los calcularemos con los datos de entrenamiento y luego veremos su rendimiento con los datos de validación. Tomar un 80% es relativamente arbitrario.
El modelo de regresión lineal con todas las variables - PRIMER INTENTO
En este primer intento, vamos a incluir todas las variables en el modelo. Dado que tenemos pocas variables, es siempre mejor empezar incluyendo todas y quitar de una en una.
modelo<-lm(Sales ~ . ,data=train)
summary(modelo)##
## Call:
## lm(formula = Sales ~ ., data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9995 -0.8527 0.0767 0.8563 3.7891
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3999017 0.3303581 13.319 <2e-16 ***
## TV 0.0561753 0.0014817 37.913 <2e-16 ***
## Radio 0.1049144 0.0090859 11.547 <2e-16 ***
## Newspaper 0.0006471 0.0061189 0.106 0.916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.606 on 156 degrees of freedom
## Multiple R-squared: 0.9138, Adjusted R-squared: 0.9121
## F-statistic: 551 on 3 and 156 DF, p-value: < 2.2e-16Claramente, con un p-valor de 0.9, sobra la variable Newspaper de nuestro modelo lineal. Los anuncios de prensa escrita no son un buen predictor. Esto era lo esperado, dados los gráficos arriba y el coeficiente de correlación en torno a 0.2 de Sales con Newspaper.
Observar que el coeficiente de determinación es muy alto, indicando un buen ajuste.
El modelo de regresión lineal con TV y Radio
Quitamos la variable Newspaper
modelo<-lm(Sales ~ TV + Radio ,data=train)
summary(modelo)##
## Call:
## lm(formula = Sales ~ TV + Radio, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.0188 -0.8520 0.0761 0.8500 3.7813
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.410791 0.312915 14.10 <2e-16 ***
## TV 0.056180 0.001476 38.05 <2e-16 ***
## Radio 0.105270 0.008415 12.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.601 on 157 degrees of freedom
## Multiple R-squared: 0.9138, Adjusted R-squared: 0.9127
## F-statistic: 831.8 on 2 and 157 DF, p-value: < 2.2e-16Ahora sí, todas las variables son significativas. Notar también que el coeficiente de determinación, de nuevo, es muy alto, indicando un muy buen ajuste al modelo.
Estos son los coeficientes del hiperplano de regresión:
modelo$coefficients## (Intercept) TV Radio
## 4.41079054 0.05617971 0.10526967Que se traduce en Sales = 4.41079054 + 0.05617971 TV + 0.10526967 Radio
La validación del modelo
Veamos cuánto difieren las predicciones de nuestro segundo modelo, con sólo las variables Radio y TV, de los datos que ya conocemos y que hemos guardado en la partición “test”.
#cálculo de predicciones
library(caret)
predictions<-predict(modelo,test)Si lo plasmamos en un gráfico:
plot(test$Sales, predictions, main="ventas vs predicciones")
abline(coef = c(0,1), col="red")
Los punto deberían caer sobre la diagonal y, sin embargo, conforme
aumenta el valor de Sales van siendo cada vez peores. Aquí tenemos un
primer vistazo a la media y la desviación típica de los errores que
hemos cometido, que parecen querer decir que el modelo no ajusta.
mean(predictions - advertising$Sales)## [1] -0.1850517sd(predictions - advertising$Sales)## [1] 7.405699Teniendo en cuenta que Sales tiene un valor máximo de 27, una desviación típica de 7.40 es bastante grande.
¿Qué está sucediendo? - Comprobando las hipótesis de la regresión lineal
Hay varias cosas que deberemos validar antes de dar el modelo por bueno. Podemos hacerlo gráficamente - al menos en principio - o usando diversos test de hipótesis.
#leemos los residuos
res <- resid(modelo)Lo primero que haremos será el gráfico con los residuos ajustados:
plot(fitted(modelo), res)
abline(h=0,col="red")
El modelo parece lineal y los residuos tienen pinta de estar
distribuídos aleatoriamente en torno a la recta y=0, pero no es tan
sencillo determinar si tiene normalidad u homocedasticidad a simple
vista. Es más, si vemos el qqplot para determinar la normalidad de los
residuos nos encontraremos con una desagradable sorpresa:
qqnorm(res)
qqline(res) 
Esto no tiene aspecto de ser normal. Debemos recurrir a test estadísticos para confirmar la validez del modelo.
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.97622, p-value = 0.007321Y ahora sí, tenemos la confirmación definitiva de que no se cumple la normalidad de los residuos.
Homocedasticidad
En el caso de la homocedasticidad para los residuos,vamos a usar el test de Breush-Pagan. La Hipótesis Nula es la homocedasticidad de los residuos y la alternativa la heterocedasticidad:
library(lmtest)
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 0.44295, df = 2, p-value = 0.8013La homocedasticidad, sin embargo, si que se está cumpliendo.
¿Es muy grave que no se cumpla la hipótesis de normalidad?
En nuestro caso, hemos visto que la dispersión de los valores es mayor de lo que cabría esperar. En general, el no cumplirse la normalidad lleva predicciones cuyo ajuste no se puede garantizar. Pero esto no descarta el modelo por completo.
De hecho, hemos conseguido un ajuste decente y hemos obtenido bastante información del tipo de relación que tiene cada una de las variables independientes con Sales, con lo que podremos construir un modelo mejor a partir de éste.
Comentarios
Publicar un comentario