0. INTRODUCCION:
Una cadena de comida rapida planea agregar un nuevo articulo a su menu. Sin embargo, todavia estan indecisos entre tres posibles campañas de marketing para promocionar el nuevo producto. Para determinar que promocion tiene el mayor efecto en las ventas, el nuevo articulo se introduce en ubicaciones en varios mercados seleccionados al azar. Se utiliza una promocion diferente en cada ubicacion y las ventas semanales del nuevo articulo se registran durante las primeras cuatro semanas. Fuente:https://www.kaggle.com/chebotinaa/fast-food-marketing-campaign-ab-test?select=WA_Marketing-Campaign.csv
1. CARGA Y PREPARACION DE DATOS
library(tidyverse)
library(knitr)
library(kableExtra)df <- read.csv("WA_Marketing-Campaign.csv")
str(df)## 'data.frame': 548 obs. of 7 variables:
## $ MarketID : int 1 1 1 1 1 1 1 1 1 1 ...
## $ MarketSize : chr "Medium" "Medium" "Medium" "Medium" ...
## $ LocationID : int 1 1 1 1 2 2 2 2 3 3 ...
## $ AgeOfStore : int 4 4 4 4 5 5 5 5 12 12 ...
## $ Promotion : int 3 3 3 3 2 2 2 2 1 1 ...
## $ week : int 1 2 3 4 1 2 3 4 1 2 ...
## $ SalesInThousands: num 33.7 35.7 29 39.2 27.8 ...kable_styling( kable(head(df)))| MarketID | MarketSize | LocationID | AgeOfStore | Promotion | week | SalesInThousands |
|---|---|---|---|---|---|---|
| 1 | Medium | 1 | 4 | 3 | 1 | 33.73 |
| 1 | Medium | 1 | 4 | 3 | 2 | 35.67 |
| 1 | Medium | 1 | 4 | 3 | 3 | 29.03 |
| 1 | Medium | 1 | 4 | 3 | 4 | 39.25 |
| 1 | Medium | 2 | 5 | 2 | 1 | 27.81 |
| 1 | Medium | 2 | 5 | 2 | 2 | 34.67 |
Seleccionamos los elementos factor y variable cuantitativa que vamos a escoger de nuestros datos para el analisis:
df1 <- df%>%
select(Promotion, SalesInThousands)
kable_styling( kable(head(df1)))| Promotion | SalesInThousands |
|---|---|
| 3 | 33.73 |
| 3 | 35.67 |
| 3 | 29.03 |
| 3 | 39.25 |
| 2 | 27.81 |
| 2 | 34.67 |
Convertimos a factor la columna de Promociones y desplegamos la tabla de contigencia con el numero de observaciones por grupo o promocion:
df1$Promotion <- as.factor(df1$Promotion)
str(df1)## 'data.frame': 548 obs. of 2 variables:
## $ Promotion : Factor w/ 3 levels "1","2","3": 3 3 3 3 2 2 2 2 1 1 ...
## $ SalesInThousands: num 33.7 35.7 29 39.2 27.8 ...table.promocion <- table(df1$Promotion)
table.promocion##
## 1 2 3
## 172 188 1882. VISUALIZACION ESTADISTICA Y GRAFICA GLOBAL DE NUESTROS DATOS
Hacemos una primera visualizacion de nuestros datos a nivel descriptivo:
library(psych)describeBy(df1$SalesInThousands, df1$Promotion)##
## Descriptive statistics by group
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 172 58.1 16.55 55.39 56.49 13.22 30.81 99.65 68.84 0.85 -0.01
## se
## X1 1.26
## ------------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 188 47.33 15.11 45.39 45.85 10.34 17.34 88.64 71.3 0.91 0.62
## se
## X1 1.1
## ------------------------------------------------------------
## group: 3
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 188 55.36 16.77 51.16 54.14 12.25 22.18 96.48 74.3 0.76 -0.21
## se
## X1 1.22Posteriormente hacemos un boxplot para tener una vision global grafica de nuestros datos:
ggplot(data=df1,aes(x =Promotion, y = SalesInThousands, color = Promotion)) +
geom_boxplot()+
theme_bw()
library(rstatix)Los valores atipicos se pueden identificar facilmente utilizando metodos de diagrama de caja, implementados en la funcion R identify_outliers()[paquete rstatix].
outliers <- df1 %>%
group_by(Promotion) %>%
identify_outliers(SalesInThousands)
kable_styling( kable(outliers))| Promotion | SalesInThousands | is.outlier | is.extreme |
|---|---|---|---|
| 1 | 94.17 | TRUE | FALSE |
| 1 | 93.71 | TRUE | FALSE |
| 1 | 96.01 | TRUE | FALSE |
| 1 | 93.03 | TRUE | FALSE |
| 1 | 97.61 | TRUE | FALSE |
| 1 | 94.43 | TRUE | FALSE |
| 1 | 91.60 | TRUE | FALSE |
| 1 | 93.86 | TRUE | FALSE |
| 1 | 99.65 | TRUE | FALSE |
| 1 | 99.12 | TRUE | FALSE |
| 1 | 93.32 | TRUE | FALSE |
| 1 | 91.29 | TRUE | FALSE |
| 2 | 87.43 | TRUE | FALSE |
| 2 | 81.79 | TRUE | FALSE |
| 2 | 88.12 | TRUE | FALSE |
| 2 | 75.29 | TRUE | FALSE |
| 2 | 88.64 | TRUE | FALSE |
| 2 | 81.37 | TRUE | FALSE |
| 2 | 82.14 | TRUE | FALSE |
| 2 | 79.64 | TRUE | FALSE |
| 2 | 73.22 | TRUE | FALSE |
| 2 | 75.88 | TRUE | FALSE |
| 2 | 78.01 | TRUE | FALSE |
| 2 | 80.17 | TRUE | FALSE |
| 2 | 82.65 | TRUE | FALSE |
| 2 | 77.39 | TRUE | FALSE |
| 2 | 80.83 | TRUE | FALSE |
| 2 | 80.75 | TRUE | FALSE |
| 2 | 82.86 | TRUE | FALSE |
| 2 | 83.40 | TRUE | FALSE |
| 2 | 75.61 | TRUE | FALSE |
| 2 | 79.53 | TRUE | FALSE |
| 2 | 74.03 | TRUE | FALSE |
| 2 | 78.53 | TRUE | FALSE |
| 2 | 76.71 | TRUE | FALSE |
| 2 | 17.34 | TRUE | FALSE |
| 3 | 89.70 | TRUE | FALSE |
| 3 | 90.30 | TRUE | FALSE |
| 3 | 89.77 | TRUE | FALSE |
| 3 | 88.91 | TRUE | FALSE |
| 3 | 94.21 | TRUE | FALSE |
| 3 | 96.48 | TRUE | FALSE |
| 3 | 91.98 | TRUE | FALSE |
| 3 | 94.89 | TRUE | FALSE |
| 3 | 93.63 | TRUE | FALSE |
| 3 | 91.61 | TRUE | FALSE |
y vemos que de los outliers de nuestra relacion no hay ninguno extremo.
3. ANOVA DE UNA VIA O FACTOR
Donde haremos comparaciones multiples para evaluar la diferencia de medias
CONTRASTE DE HIPOTESIS:
H0: la media de todas las muestras son iguales
H1: al menos una de las medias es diferente del resto
p>0.05 aceptamos la hipotesis nula H0
p<0.05 aceptamos la hipotesis alternativa H1El Anova requiere el cumplimiento los siguientes supuestos:
-Las poblaciones (distribuciones de probabilidad de la variable dependiente correspondiente a cada factor) son normales.
-Las K muestras sobre las que se aplican los tratamientos son independientes.
-Las poblaciones tienen todas igual varianza (homocedasticidad).
Comprobamos ahora si se pueden cumplir los supuestos para realizar el ANOVA:
GRAFICAS PARA VER SI SIGUE UNA DISTRIBUCION NORMAL
par(mfrow = c(2,2))
qqnorm(df1[df1$Promotion == "1","SalesInThousands"], main = "1")
qqline(df1[df1$Promotion == "1","SalesInThousands"])
qqnorm(df1[df1$Promotion == "2","SalesInThousands"], main = "2")
qqline(df1[df1$Promotion == "2","SalesInThousands"])
qqnorm(df1[df1$Promotion == "3","SalesInThousands"], main = "3")
qqline(df1[df1$Promotion == "3","SalesInThousands"])
Se puede apreciar en el grafico que nuestros datos se desvian de la linea normal en un tramo de los mismos en los distintos grupos, pero debemos comprobarlo mas feacientemente con un test de normalidad.
Los datos son independientes, procedemos pues ha hacer la prueba de normalidad de Kolmogorov-Smirnov con la correccion de Lilliefors,ya que tenemos mas de 50 observaciones en algunos grupos,en este caso todos los niveles,si fueran menores de 50 por grupo o nivel se emplearia el test de Shapiro-Wilk.
3.1 PRUEBA DE NORMALIDAD DE KOLMOGOROV-SMIRNOV (CON CORRECCION LILLIEFORS)
CONTRASTE DE HIPOTESIS:
H0: los datos proceden de una distribucion normal
H1: los datos no proceden de una distribucion normal
p>0.05 aceptamos la hipotesis nula H0
p<0.05 aceptamos la hipotesis alternativa H1library(nortest)by(data = df1,INDICES = df1$Promotion,FUN = function(x){ lillie.test(x$SalesInThousands)})## df1$Promotion: 1
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$SalesInThousands
## D = 0.12268, p-value = 1.235e-06
##
## ------------------------------------------------------------
## df1$Promotion: 2
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$SalesInThousands
## D = 0.14253, p-value = 5.672e-10
##
## ------------------------------------------------------------
## df1$Promotion: 3
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$SalesInThousands
## D = 0.14589, p-value = 1.778e-10Como se puede ver en los 3 grupos de campañas promocionales el pvalor<0.05, lo cual hace que se rechace la hipotesis de que sean normales estas distribuciones,al no ser normal inclumple uno de los supuestos de ANOVA.
Procedemos a probar ahora la homocedasticidad(igualdad de varianzas):
3.2 PRUEBA DE HOMOCEDASTICIDAD DE BARTLETT
CONTRASTE DE HIPOTESIS:
H0: las muestras presentan varianzas iguales
H1: las muestras presentan varianzas distintas
p>0.05 aceptamos la hipotesis nula H0
p<0.05 aceptamos la hipotesis alternativa H1bartlett.test(df1$SalesInThousands~df1$Promotion)##
## Bartlett test of homogeneity of variances
##
## data: df1$SalesInThousands by df1$Promotion
## Bartlett's K-squared = 2.3332, df = 2, p-value = 0.3114El resultado nos da un pvalor>0.05, por tanto nos da que no hay evidencias en contra de la homogeneidad de varianzas.
Una vez que hemos visto que no se cumplen integramente los supuestos ANOVA, la prueba alternativa no parametrica es la de el test de Kruskal-Wallis, siempre que cumpla una serie de condiciones:
-No es necesario que las muestras que se comparan provengan de una distribucion normal.
-Homocedasticidad: dado que la hipotesis nula asume que todos los grupos pertenecen a una misma poblacion y que por lo tanto tienen las mismas medianas, es necesario que todos los grupos tengan la misma varianza.se puede comprobar con el test de Levene o Bartlett.
-Misma distribucion para todos los grupos: la distribucion de los grupos no tiene que ser normal, pero ha de ser igual en todos(por ejemplo que todos muestren asimetria hacia la izquierda).
Examinamos pues estas condiciones:
-Efectivamente no tiene por que provenir de una distribucion normal, como es nuestro caso.
-Vamos a ver si hay semejanza entre la distribucion de grupos.
ggplot(data = df1, aes(x=SalesInThousands, colour = Promotion))+
geom_histogram() +
theme_bw() +
facet_grid(.~ Promotion )+
theme(legend.position = "none")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
-Las 3 distribuciones parecen iguales o muy parecidas como puede verse en los graficos respectivos.
-y por ultimo, como ya hemos visto con el test de Bartlett, no hay evidencias en contra de la homogeneidad de varianzas (homocedasticidad).
Por tanto como cumple las condiciones, vamos a hacer la prueba en sustitucion de la ANOVA,La prueba no parametrica de Kruskal-Wallis.
4. TEST DE KRUSKAL-WALLIS
A diferencia del ANOVA en el que se comparan medias, el test de Kruskal-Wallis contrasta si las diferentes muestras estan equidistribuidas y que por lo tanto pertenecen a una misma distribucion (poblacion). Bajo ciertas simplificaciones puede considerarse que el test de Kruskal-Wallis compara las medianas.
CONTRASTE DE HIPOTESIS:
H0: Med1= Med2=…=Medk
H1: Medi ≠ Medj al menos para un par (i,j)
H0: Las muestras provienen de poblaciones identicas
H1: Las muestras provienen de poblaciones diferentes
p>0.05 aceptamos la hipotesis nula H0
p<0.05 aceptamos la hipotesis alternativa H1kruskal.test(SalesInThousands~Promotion, data=df1)##
## Kruskal-Wallis rank sum test
##
## data: SalesInThousands by Promotion
## Kruskal-Wallis chi-squared = 53.295, df = 2, p-value = 2.674e-12En este caso el pvalor<0.05, por tanto hay evidencia de que al menos hay diferencias de al menos 2 grupos.
4.1 COMPARACIONES POST HOC PARA SABER QUE PARES DE GRUPOS DIFIEREN
De entre los 2 metodos de comparacion mas empleados: Test de Mann-Whitney y Tukey´s range test,y de entre los diferentes metodos de correccion de significancia entre los que destacan “Bonferroni” y “Holm”, vamos a elegir el de “Holm”.
4.1.1 TEST DE MANN-WHITNEY
CONTRASTE DE HIPOTESIS:
H0:No hay diferencias entre grupos
H1:Si hay diferencias entre grupos
p>0.05 aceptamos la hipotesis nula H0
p<0.05 aceptamos la hipotesis alternativa H1Condiciones de cumplimiento para realizar el test:
-Los datos tienen que ser independientes.
-Los datos tienen que ser ordinales o bien se tienen que poder ordenar de menor a mayor.
-No es necesario asumir que las muestras se distribuyen de forma normal o que proceden de poblaciones normales. Pero, para que el test compare medianas, ambas han de tener el mismo tipo de distribucion (varianza, asimetria).
-Igualdad de varianza entre grupos (homocedasticidad).
Como se ve cumple con todas las condiciones, asi que continuamos:
pairwise.wilcox.test(df1$SalesInThousands,df1$Promotion, p.adjust = "holm")##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: df1$SalesInThousands and df1$Promotion
##
## 1 2
## 2 1.8e-11 -
## 3 0.035 2.4e-07
##
## P value adjustment method: holmComo se puede ver en todas las comparaciones por pares, se ve que el pvalor<0.05, eso quiere decir que se rechaza la hipotesis nula y se demuestra que hay evidencias de que las medianas son diferentes entre los 3 grupos sometidos a estudio.
Como se puede ver en los descriptivos a continuacion y teniendo en cuenta en el resultado del test anterior que las medianas son significativamente diferentes de cada promocion entre si, tenemos:
results <- aggregate(SalesInThousands ~ Promotion, data = df1, FUN = median)
kable_styling( kable(results))| Promotion | SalesInThousands |
|---|---|
| 1 | 55.385 |
| 2 | 45.385 |
| 3 | 51.165 |
Tambien podemos visualizar estas diferencias en el grafico de barras donde se representa la mediana de los datos y la prevalencia de la promocion 1 con respecto a las demas.
ggplot(data=results,aes(Promotion, SalesInThousands, fill=Promotion)) +
geom_col()+
geom_text(position="identity",vjust=0.02,aes(Promotion, SalesInThousands, label=sprintf("%0.2f", round(SalesInThousands, digits = 2))))
Por tanto promocion1>promocion3>promocion2. La promocion1 tiene mayor efecto sobre las ventas siendo los grupos entre si diferentes.