En esta práctica voy a trabajar con el dataset marketing del paquete datarium. Esta base de datos está formada por 200 observaciones con 4 variables. Estudia la inversión en publicidad en tres plataformas diferentes (youtube, facebook y newspaper) y su relación con las ventas.
En primer lugar compruebo el tipo de datos y la existencia de datos perdidos:
Todas las variables son contínuas, y no tenemos datos perdidos. El siguiente panel muestra la relación entre las diferentes variables, así como la distribución de cada una de ellas:
Vemos que existe una correlación significativa entre las ventas y la inversión en publicidad en los diferentes medios. Además también destaca la correlación (r=0.35) entre la inversión en newspaper y facebook. Podemos observar también que los datos de las diferentes variables no siguen una distribución normal, excepto la variable sales. El siguiente gráfico muestra mejor la fuerza de la correlación entre las diferentes variables:
Divido aleatoriamente el dataset en un grupo de entrenamiento (80%) y otro grupo test (20%). En primer lugar realizo un test de significación global, comparando el modelo completo vs el modelo nulo. Si el test es no significativo no hará falta seguir con el análisis:
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 159 | 5680.4464 | NA | NA | NA | NA |
| 156 | 678.8336 | 3 | 5001.613 | 383.1334 | 0 |
Así pues confirmamos que existen diferencias entre el modelo completo y el modelo nulo. La siguiente cuestión a responder es cúal es el mejor modelo .
El summary del modelo completo lo muestro a continuación. Comprobamos que la variable newspaper es no significativa.
##
## Call:
## lm(formula = sales ~ ., data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.3187 -1.0570 0.3526 1.4679 3.3944
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.527417 0.442016 7.980 2.95e-13 ***
## youtube 0.046232 0.001633 28.316 < 2e-16 ***
## facebook 0.183196 0.009888 18.526 < 2e-16 ***
## newspaper -0.003166 0.006638 -0.477 0.634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.086 on 156 degrees of freedom
## Multiple R-squared: 0.8805, Adjusted R-squared: 0.8782
## F-statistic: 383.1 on 3 and 156 DF, p-value: < 2.2e-16
El contraste entre el modelo completo y el modelo sin newspaper lo calculo con el siguiente anova:
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 157 | 679.8238 | NA | NA | NA | NA |
| 156 | 678.8336 | 1 | 0.9901896 | 0.2275515 | 0.6340129 |
Puesto que el modelo reducido presenta un RSE ligeramente menor (2.018 vs 2.023) con un R2 igual (0.89), me quedo con el modelo reducido por ser el más sencillo.
1- Normalidad de los residuos. Los siguientes gráficos muestran cómo los residuales se alejan de la distribución normal. El test de Shapiro-Wilks confirma este hallazgo.
##
## Shapiro-Wilk normality test
##
## data: residuals(model.1)
## W = 0.90939, p-value = 2.078e-08
##
## studentized Breusch-Pagan test
##
## data: model.1
## BP = 2.9425, df = 2, p-value = 0.2296
| sales | fb | yt | ti | res | |
|---|---|---|---|---|---|
| 131 | 1.92 | 47.52 | 0.84 | 5.419857 | -10.214159 |
| 6 | 8.64 | 58.68 | 10.44 | 2.999541 | -5.964111 |
| 179 | 14.16 | 2.76 | 332.04 | 2.556049 | -5.140642 |
En cuanto al impacto sobre los coeficientes (dfbetas), ninguna de las observaciones supera el umbral de 28.2. El valor máximo dfbeta lo muestro a continuación:
## [1] 0.1459795
El impacto sobre el ajuste lo mido mediante el parámetro dffit. Utilizo un punto de corte de 0.244. A continuación muestro las observaciones que lo superan:
## 131 6 36 179 127 26 76 132
## 0.9780048 0.6159521 0.4790720 0.4719195 0.3814500 0.3728543 0.3110953 0.3059703
## 79 3 103 159 118 148 77 189
## 0.3039976 0.3030858 0.2779556 0.2754132 0.2675416 0.2664235 0.2639602 0.2594268
## 57 184
## 0.2567373 0.2543682
No obstante, a pesar de que existan observaciones influyentes las mantengo en el dataset puesto que considero que son datos que han sido tomados correctamente y forman parte de la variabilidad del modelo.
Puesto que el criterio de normalidad no se cumple, utilizaré una regresión robusta, menos sensible a la influencia de datos atípicos. El resumen del modelo lo muestro a continuación:
##
## Call: rlm(formula = sales ~ facebook + youtube, data = train)
## Residuals:
## Min 1Q Median 3Q Max
## -10.8057 -1.2581 0.2507 1.2136 3.2211
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 3.7857 0.3924 9.6486
## facebook 0.1873 0.0086 21.6979
## youtube 0.0446 0.0015 29.5325
##
## Residual standard error: 1.857 on 157 degrees of freedom
Si lo comparamos con el modelo reducido calculado de forma ordinaria (OLS), vemos que no hay mucha diferencia en el valor de los coeficientes de regresión en este caso:
##
## Call:
## lm(formula = sales ~ . - newspaper, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.2142 -0.9847 0.3927 1.5088 3.4060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.465865 0.421721 8.218 7.23e-14 ***
## youtube 0.046180 0.001625 28.418 < 2e-16 ***
## facebook 0.181597 0.009280 19.568 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.081 on 157 degrees of freedom
## Multiple R-squared: 0.8803, Adjusted R-squared: 0.8788
## F-statistic: 577.4 on 2 and 157 DF, p-value: < 2.2e-16
El coeficiente R2 para el grupo de entrenamiento del modelo robusto sería:
## [1] 0.879443
Finalmente valido el modelo en el grupo test. El error de validación en el grupo test (RMSE) lo calculo a continuación:
## [1] 1.218909
De nuevo vemos que es muy parecido al error de validación mediante regresión ordinaria:
## [1] 1.216325
Aplico el modelo al dataset global, y así obtengo los coeficientes definitivos de la recta de regresión. En el siguiente gráfico muestro la similitud que para este caso tienen las dos rectas, regresión robusta y regresión OLS:
## Warning in abline(model.1, col = "red", lwd = 2): only using the first two of 3
## regression coefficients
## Warning in abline(rmodel, col = "blue", lwd = 2): only using the first two of 3
## regression coefficients