Librerías que hay que instalar para seguir esta práctica:

# Esta es la librería para insertar un video de youtube
library(vembedr)

Un tutorial para instalar R-Studio

embed_url("https://youtu.be/vnY4QV-Pq6c?si=L_-nzenuKciePWWK")

Un tutorial para instalar paquetes

embed_url("https://youtu.be/zqCQxWcWYYs?si=kdjDsqss9AmXrFr1")

Iniciar el trabajo con R Markdown

embed_url("https://youtu.be/rCP7OAyPMyY?si=5qiS12UpwpzHv76W")

Ahora cargamos las librerías necesarias para esta práctica:

library(rmarkdown)
library(DataExplorer)
library(Hmisc)
library(stats)
library(dplyr)
library(kableExtra)
library(knitr)
library(ggplot2)
library(car)
library(psych)
library(GGally)
library(lsr)
library(ggpubr)
library(phia)
library(arules)
library(agricolae)

Se carga la base de datos adjunta en formato .csv

Los datos que se utilizarán para esta actividad corresponde a las ventas de asientos de coches infantiles a 400 tiendas distintas en US. Las variables son: • Sales (Ventas unitarias, en miles, en cada ubicación) • CompPrice (Precio cobrado por el competidor en cada ubicación) • Income (Nivel de ingresos comunitarios, en miles de dólares) • Advertising (Presupuesto de publicidad local de la empresa en cada ubicación, en miles de dólares) • Population (Tamaño de la población en la región, en miles) • Price (Precio para asientos de coches en cada sitio) • ShelveLoc (Un factor con niveles Bad, Good y Medium que indica la calidad de la ubicación de los estantes de los asientos del coche de cada sitio) • Age (Edad media de la población local) • Education (Nivel educativo en cada sitio) • Urban (Un factor con los niveles Yes y No para indicar si la tienda se encuentra en una ubicación urbana o rural) • US (Un factor con los niveles Yes y No para indicar si la tienda se encuentra en EE.UU. o no) Los datos del estudio están en el archivo ChildCarSeats1.csv

Cualquier archivo excel se puede convertir mediante “guardar como” en formato .csv. Para cargar un archivo en formato SPSS, seguir este tutorial.

df<- read.csv('ChildCarSeats1.csv')

1 Estadística descriptiva y visualización

1.1 Tipo de datos

utilizamos el método ‘summary’

summary(df)
##      Sales          CompPrice       Income        Advertising    
##  Min.   : 0.000   Min.   : 77   Min.   : 21.00   Min.   : 0.000  
##  1st Qu.: 5.390   1st Qu.:115   1st Qu.: 42.75   1st Qu.: 0.000  
##  Median : 7.435   Median :125   Median : 69.00   Median : 5.000  
##  Mean   : 7.410   Mean   :125   Mean   : 68.66   Mean   : 6.635  
##  3rd Qu.: 9.160   3rd Qu.:135   3rd Qu.: 91.00   3rd Qu.:12.000  
##  Max.   :16.270   Max.   :175   Max.   :120.00   Max.   :29.000  
##    Population        Price        ShelveLoc              Age       
##  Min.   : 10.0   Min.   : 24.0   Length:400         Min.   :25.00  
##  1st Qu.:139.0   1st Qu.:100.0   Class :character   1st Qu.:39.75  
##  Median :272.0   Median :117.0   Mode  :character   Median :54.50  
##  Mean   :264.8   Mean   :115.8                      Mean   :53.32  
##  3rd Qu.:398.5   3rd Qu.:131.0                      3rd Qu.:66.00  
##  Max.   :509.0   Max.   :191.0                      Max.   :80.00  
##    Education       Urban                US           
##  Min.   :10.0   Length:400         Length:400        
##  1st Qu.:12.0   Class :character   Class :character  
##  Median :14.0   Mode  :character   Mode  :character  
##  Mean   :13.9                                        
##  3rd Qu.:16.0                                        
##  Max.   :18.0

Donde nos indica que existen 11 variables: Ya vemos que hay 3 que son categóricas: “ShelveLoc”, “Urbarn” y “US”. El resto son variables continuas.

str(df)
## 'data.frame':    400 obs. of  11 variables:
##  $ Sales      : num  9.5 11.22 10.06 7.4 4.15 ...
##  $ CompPrice  : int  138 111 113 117 141 124 115 136 132 132 ...
##  $ Income     : int  73 48 35 100 64 113 105 81 110 113 ...
##  $ Advertising: int  11 16 10 4 3 13 0 15 0 0 ...
##  $ Population : int  276 260 269 466 340 501 45 425 108 131 ...
##  $ Price      : int  120 83 80 97 128 72 108 120 124 124 ...
##  $ ShelveLoc  : chr  "Bad" "Good" "Medium" "Medium" ...
##  $ Age        : int  42 65 59 55 38 78 71 67 76 76 ...
##  $ Education  : int  17 10 12 14 13 16 15 10 10 17 ...
##  $ Urban      : chr  "Yes" "Yes" "Yes" "Yes" ...
##  $ US         : chr  "Yes" "Yes" "Yes" "Yes" ...

También podemos ayudar de a la función “describe” del este Hmisc.

describe(df)
##             vars   n   mean     sd median trimmed    mad min    max  range
## Sales          1 400   7.41   2.73   7.44    7.36   2.79   0  16.27  16.27
## CompPrice      2 400 124.97  15.33 125.00  125.04  14.83  77 175.00  98.00
## Income         3 400  68.66  27.99  69.00   68.26  35.58  21 120.00  99.00
## Advertising    4 400   6.64   6.65   5.00    5.89   7.41   0  29.00  29.00
## Population     5 400 264.84 147.38 272.00  265.56 191.26  10 509.00 499.00
## Price          6 400 115.80  23.68 117.00  115.92  22.24  24 191.00 167.00
## ShelveLoc*     7 400   2.31   0.83   3.00    2.38   0.00   1   3.00   2.00
## Age            8 400  53.32  16.20  54.50   53.48  20.02  25  80.00  55.00
## Education      9 400  13.90   2.62  14.00   13.88   2.97  10  18.00   8.00
## Urban*        10 400   1.71   0.46   2.00    1.76   0.00   1   2.00   1.00
## US*           11 400   1.65   0.48   2.00    1.68   0.00   1   2.00   1.00
##              skew kurtosis   se
## Sales        0.14    -0.06 0.14
## CompPrice   -0.04     0.01 0.77
## Income       0.05    -1.10 1.40
## Advertising  0.63    -0.57 0.33
## Population  -0.05    -1.21 7.37
## Price       -0.12     0.41 1.18
## ShelveLoc*  -0.62    -1.28 0.04
## Age         -0.08    -1.14 0.81
## Education    0.04    -1.31 0.13
## Urban*      -0.90    -1.20 0.02
## US*         -0.60    -1.64 0.02

Nos ofrece más información con proporciones. Tenemos pues 11 variables, 3 de tipo Factor, 8 variables continuas. Análisis de los valores faltantes:

plot_missing(df)

NO existen “missing values”.

Definimos las variables categóricas

df$ShelveLoc = as.factor(df$ShelveLoc)
df$US = as.factor(df$US)
df$Urban = as.factor(df$Urban)

2 Estadística descriptiva i visualitzación

Seleccionamos las variables continuas y creamos una dataset aparte. Nos fijamos que las columnas 7,10 y 11 (en el orden) las quitamos

df_c <- df[-c(7,10,11)]

Las columnas iniciales, que son 11

colnames(df)
##  [1] "Sales"       "CompPrice"   "Income"      "Advertising" "Population" 
##  [6] "Price"       "ShelveLoc"   "Age"         "Education"   "Urban"      
## [11] "US"

Si quitamos las columnas 7,10,11

colnames(df_c) # Muestra el nombre de las columnas
## [1] "Sales"       "CompPrice"   "Income"      "Advertising" "Population" 
## [6] "Price"       "Age"         "Education"

Vamos a crear una tabla que nos muestre 4 variables estadísticas:

  • La media aritmética
  • La mediana
  • La desviacion standar
  • La amplitud intercuartílica (IQR en inglés)

Voy comentando el código

media <- vector()   #Genero los vectores para almacenar los resultados 
mediana <- vector() # Uno por cada variable que analizo.
sd <- vector()
IQR <- vector()

for (i in 1:dim(df_c)[2]){          # Genero un loop
  media[i]= mean(df_c[,i])          # Recorro y guardo cada estadístico de cada variable contínua
  mediana[i]= median(df_c[,i])
  sd[i] = sd(df_c[,i])
  IQR[i]= IQR(df_c[,i])
}

Qué contiene cada vector? Los resultados de cada variable contínua (las 7 que hemos seleccionado)

media
## [1]   7.409975 124.975000  68.657500   6.635000 264.840000 115.795000  53.322500
## [8]  13.900000

Otra forma más sencilla de conseguir estos datos

describe(df_c)
##             vars   n   mean     sd median trimmed    mad min    max  range
## Sales          1 400   7.41   2.73   7.44    7.36   2.79   0  16.27  16.27
## CompPrice      2 400 124.97  15.33 125.00  125.04  14.83  77 175.00  98.00
## Income         3 400  68.66  27.99  69.00   68.26  35.58  21 120.00  99.00
## Advertising    4 400   6.64   6.65   5.00    5.89   7.41   0  29.00  29.00
## Population     5 400 264.84 147.38 272.00  265.56 191.26  10 509.00 499.00
## Price          6 400 115.80  23.68 117.00  115.92  22.24  24 191.00 167.00
## Age            7 400  53.32  16.20  54.50   53.48  20.02  25  80.00  55.00
## Education      8 400  13.90   2.62  14.00   13.88   2.97  10  18.00   8.00
##              skew kurtosis   se
## Sales        0.14    -0.06 0.14
## CompPrice   -0.04     0.01 0.77
## Income       0.05    -1.10 1.40
## Advertising  0.63    -0.57 0.33
## Population  -0.05    -1.21 7.37
## Price       -0.12     0.41 1.18
## Age         -0.08    -1.14 0.81
## Education    0.04    -1.31 0.13

De hecho, si selecciono sólo la columna “mean” que significa media aritmética:

d = describe(df_c) # Convierto estos datos en una tabla (base de datos o dataframe)
d$mean             # ahora sólo muestro los datos de la media aritmética, Utilizo $ para seleccionar las columnas
## [1]   7.409975 124.975000  68.657500   6.635000 264.840000 115.795000  53.322500
## [8]  13.900000
tabla_c<-list(colnames(df_c),media, mediana, sd, IQR)      # Las guardo en una tabla
tabla_c<-as.data.frame(tabla_c, col.names = c('variables', 'media','mediana','variancia','IQR')) # les doy nombre a cada columna

2.1 Resumen datos cuantitativos

2.1.1 Estudio de la variabilidad

Una vez con la tabla anterior, podemos realizar un estudio sobre la variabilidad de los datos. Podemos completar la tabla con el coeficiente de variación (ratientre desviación típica y media).

tabla_c$cv <- tabla_c$variancia/tabla_c$media
tabla_c %>%
  kable() %>%
  kable_styling()
variables media mediana variancia IQR cv
Sales 7.409975 7.435 2.732580 3.77 0.3687705
CompPrice 124.975000 125.000 15.334511 20.00 0.1227006
Income 68.657500 69.000 27.986037 48.25 0.4076181
Advertising 6.635000 5.000 6.650364 12.00 1.0023156
Population 264.840000 272.000 147.376436 259.50 0.5564735
Price 115.795000 117.000 23.676664 31.00 0.2044705
Age 53.322500 54.500 16.200297 26.25 0.3038173
Education 13.900000 14.000 2.620528 4.00 0.1885272

2.1.2 Comentario

Nos fijamos en el valor cv = Coeficiente de Pearson

  • Las variables que cuentan con mayor variabilidad son: Advertising - con un ratio > 1.00 Population - con un ratio de 0.55 Income - con un ratio de 0.40

la información que nos ofrece la tabla es que existe una gran diferencia entre territorios sobre la publicidad. Sabemos por las tablas anteriores que existen territorios que no han tenido gasto en publicidad y otros que han dedicado cantidades relativamente altas.

También hay que tener en cuenta que tanto la población como la renta disponible de la zona tienen una importante “anchura”, o valores altos de IQR.

Esto nos informa de que probablemente estas variables pueden tener protagonismo respecto a posibles análisis posterior. Es más probable que una variable con gran variabilidad pueda tener un elemento discriminante mayor que otros que lo tengan menor, y esto nos puede servir para predecir variables dependientes.

  • Las variables con menor variabilidad son: Precio de la competencia - con un ratio de 0.12 Nivel Educativo - con un ratio de 0.19 Price - con un ratio de 0.20

Estas variables con poca variabilidad interterritorial probablemente nos ofrecerá poca capacidad discriminante y menos información para predecir variables dependientes.

Agregar que la relación entre IQR y la distribución nos aportaría una medida similar. Un IQR de 3.77 no nos aporta información para evaluar, sino la referenciamos con otra variable, como podría ser la distribución acumulada. Nos ofrecería una información de cuanto mayor es la anchura de la distancia intercuartílica. Se espera que el ancho Q1-Q3 represente el 50% de la distribución acumulada de la variable. A mayor valor del coeficiente de variabilidad, mayor valor respecto a la media de IQR.

2.2 Diagramas de caja

2.2.1 ShelveLoc vs Salas

Dibujamos un boxplot con las tres cajas correspondientes a los tres niveles de la variable categórica ShelveLoc respecto a sus ventas. El gráfico nos monstrará la distribución de las ventas por clase, además he introducido un punto con la media para evaluar una posible kurtosis en sus distribuciones.

ggplot(df, aes(x=ShelveLoc, y=Sales))+ geom_boxplot()+                      # utilizo la libreria ggplot para dibujar gráficos          
  stat_summary(fun =mean, geom ="point", shape =19,size =2,color ="brown")  # Añado el punto de la media para evaluar la kurtosis, si el punto no está junto a la línea que divide la caja.

Claramente la clase “Good” supera con ventas al resto, y ya marcamos la clase “medium” con outliers values. Estos valores “fuera de rango” se muestran mediante puntos fuera de la caja, vemos en “medium” los puntos en la parte inferior.

2.2.2 Urban vs Sales

Repetimos el análisis

ggplot(df, aes(x=Urban, y=Sales))+ geom_boxplot()+
  stat_summary(fun =mean, geom ="point", shape =19,size =2,color ="brown")

En cambio en la variable “Urban” no hay diferencias según la gráfica, y anotar los outliers values en la clase “Yes”.

2.2.3 US vs Sales

ggplot(df, aes(x=US, y=Sales))+ geom_boxplot()+
  stat_summary(fun =mean, geom ="point", shape =19,size =2,color ="brown")

Aquí parece que sí hay diferencia entre que la tienda sea en US y las ventas. Podríamos decir que las tiendas en US tienden a vender más que el resto. Los outliers values están en la clase “NO”.

2.3 Representación de las variables cualitativas

Dibujaremos un gráfico de barras las tres variables cualitativas, que son ShelveLoc, Urban y US.

ggplot(df, aes(ShelveLoc))+ geom_bar(fill = "brown")+ theme_gray()  # nos fijamos en el método "geom_bar" para indicar que nos interesa gráfico de barras

Nos da una primera aproximación entre las ventas y el número de clase. Good tenía unas ventas medias muy superior a “Medium”. Lo anotamos. Ahora estamos interesados en saber la proporción “YES” y “NO” de las variables categóricas “Urban” y “US”

ggplot(df, aes(Urban))+ geom_bar(fill = "brown")+ theme_gray()

Mayoritariamente clase “Yes” pero según el gráfico box-plot anterior no existían diferencias respecto a las ventas.

ggplot(df, aes(US))+ geom_bar(fill = "brown")+ theme_gray()

Mayoría en tiendas en US.

3 Estadística inferencial

3.1 Intervalos de confianza de la variable Price

Creamos una función para calcular el intervalo de confianza.

calcIC <- function(m,n,s,z){                          # Entro las siguientes variables m, n, s, z
  ui <- m- z * (s/sqrt(n))
  us <- m+ z * (s/sqrt(n))
  R <- list(ui,us)
  cat("Intervalo de Confianza es (", ui,",",us,")\n")
  return(R)                                           # Me devuelve el intervalo de confianza
  }

donde la variable a estudiar:

m -> es la media.

n -> número de observaciones.

s -> desviación típica.

z -> estadístico que calcularemos bilateral con un 95% de confianza.

Para entender un aparte importante, un tutorial sencillo: https://www.statology.org/working-with-the-student-t-distribution-in-r-dt-qt-pt-rt/

La función qt devuelve el valor de la función de densidad acumulada inversa (cdf) de la distribución t de Student dada una cierta variable aleatoria x y grados de libertad df. Pero si lo analizamos en las 2 colas a un 95% de confianza, pues requerimos 0.975 (que representa un 2.5% por cola).

m <- mean(df$Price)     # Calculamnos la media de la variable Price. Recordemos que tenemos muchas maneras de calcularla, aquí mediante mean()
n <- length(df$Price)   # Número de registros mediante lenght()
s <- sd(df$Price)       # Valor de desviación típica mediante sd()
z <- qt(0.975, df = 400)     #

Ahora podemos calcular el Intervalo de Confianza de la variable Price con un 95% de confianza, de la variable Precio. Esto significa que el 95% de probabilidad es que se encuentre entre los valores indicados:

calcIC(m,n,s,z)
## Intervalo de Confianza es ( 113.4677 , 118.1223 )
## [[1]]
## [1] 113.4677
## 
## [[2]]
## [1] 118.1223

3.2 Test de comparación de dos medias.

Estudiar si el hecho de que una tienda esté fuera de US representa una diferencia de ventas significativa.

Ante todo estudiamos la variables de las ventas.

3.2.1 ESTUDIO DE NORMALIDAD

Estudiaremos con el test de Kolmogorov-Smirnov Para más info, ver: https://www.statology.org/kolmogorov-smirnov-test-r/

ks.test(x=df$Sales, mean(df$Sales), sd(df$Sales))
## Warning in ks.test.default(x = df$Sales, mean(df$Sales), sd(df$Sales)):
## Parameter(s) ignored
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  df$Sales and mean(df$Sales)
## D = 0.5075, p-value = 0.9875
## alternative hypothesis: two-sided

Con una p-value tan alta 0.95 aceptamos la hipótesis nula que la variable “Sales” procede de una distribución normal.

Con un histograma deberíamos representarlo mejor:

ggplot(data = df, aes(x = Sales)) +                           # Utilizo el paquete ggplot
  geom_histogram(aes(y = ..density.., fill = ..count..)) +    # Creo un histograma
  scale_fill_gradient(low = "yellow", high = "blue") +        # Utilizo un gradiente por colores.
  stat_function(fun = dnorm, colour = "firebrick",            # Añado la curva normal
                args = list(mean = mean(df$Sales),
                            sd = sd(df$Sales))) +
  ggtitle("Histograma que incluye la curva normal teòrica") +  # Añado un título.
  theme_bw()
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Con este histograma confirma una p-value tan alta. La distribución se ajusta mucho a la curva teórica.

Dibujamos los cuartiles teóricos: Más info ver en: https://www.statology.org/q-q-plot-r/

qqnorm(df$Sales, pch = 17, col = "brown", main="Quantiles teóricos")
qqline(df$Sales)

Prácticamente sigue la diagonal. Podemos afirmar sin duda que la variable “Sales” tiene una distribución normal.

3.2.2 Escribir la hipótesis nula y alternativa

La hipótesis nula es que no existen diferencias significativas entre las dos medias.

\(H_{o} : \mu_{1}= \mu_{2}\)

\(H_{1} : \mu_{1} \ne \mu_{2}\)\((bilateral)\)

\(\alpha = 0.95\)

3.2.3 Método del contraste entre dos medias

Lo que se quiere comparar de la variable continua “Sales” entre dos grupos, que serían la clase US-“SI” y US-“NO”. Estos dos grupos podrían considerarse poblaciones independientes, y que la variable a estudiar proviene de una distribución normal tal y como hemos acreditado.

Lo que queremos estudiar es si la diferencia observada entre las medias de ambos grupos es significativa. Y podríamos trabajar con el estadístico t-student para que cumple: - Las observaciones son independientes. - Las poblaciones se distribuyen de forma normal. - Homocedasticidad, que deberemos estudiar.

Para evaluar la homocedasticidad de ambas muestras, utilizaremos el Test de Levene. Más info: https://statologos.com/prueba-levenes-r/

leveneTest(y = df$Sales, group = df$US, center = 'median')
## Levene's Test for Homogeneity of Variance (center = "median")
##        Df F value    Pr(>F)    
## group   1  11.787 0.0006592 ***
##       398                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test de Levene sí encuentra diferencias significativas entre las variancias de ambos grupos.

Probamos pues el Test de Bartlett, ya que el número de observaciones por grupo es distinto:

UsYes <- df[df$US == 'Yes', 'Sales']
UsNo <- df[df$US == 'No' , 'Sales']
bartlett.test(list(UsYes, UsNo))
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(UsYes, UsNo)
## Bartlett's K-squared = 11.184, df = 1, p-value = 0.0008253

Aún así, sí encuentra diferencias en su homocedasticidad entre los grupos. Existe heterocedasticidad entre las muestras. Viola una de las condiciones para utilizar el t.test.

Entonces utilizaremos lo etadístico:

\(Z = \frac{\bar{\mathcal X_{1}}-\bar{\mathcal X_{2}}}{\sqrt{d_{1}+d_{2}}}\)

Calculamos \(d=\frac{s}{\sqrt{n}}\)

3.2.4 Contraste, valor crítico y p-value con 95% índice de confianza.

Efectuamos los cálculos Creamos valores de media aritmética “mean” de cada variable

m_Y <- mean(UsYes)
m_N <- mean(UsNo)
sd_Y <- sd(UsYes)
sd_N <- sd(UsNo)
n_Y <- length(UsYes)
n_N <- length(UsNo)
d_Y <- sd_Y^2/n_Y
d_N <- sd_N^2/n_N

Ahora ya podemos calcular el estadístico Z

Z <- (m_Y-m_N)/sqrt(d_Y + d_N)
Z
## [1] 4.970486
pnorm(Z)
## [1] 0.9999997

Ambos valores nos indican que debemos rechazar la hipotesis nula y afirmar que sí que existe una diferencia entre las poblaciones, en otras palabras, el hecho de tener una tienda en US influye en las ventas.

Lo dibujamos:

mu1 <- df %>%
  group_by(US) %>%
  summarise(grp.mean=mean(Sales))
ggplot(df, aes(x=Sales), main='Distribucions')+
  geom_density(aes(color =US)) +
  geom_vline(data = mu1, aes(xintercept =grp.mean,color =US),
  linetype="dashed") + theme_minimal()

3.3 Contraste no paramétrico

Como ya hemos visto antes, no pasar el test de homocedasticidad viola una de las condiciones para aplicar las test de Student.

La alternativa no paramétrica sería el test de Mann-Shitney-Wilcoxon.

wilcox.test(UsYes, UsNo)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  UsYes and UsNo
## W = 23085, p-value = 1.651e-05
## alternative hypothesis: true location shift is not equal to 0

Es inferior a p-value 0.05, por lo que confirmamos que sí es significativo. Lo podemos comparar con el resultado de t student, que deberíamos utilizar si se confirmara la no existencia de heterocedasticidad, y confirmaríamos también.

t.test(UsYes, UsNo)
## 
##  Welch Two Sample t-test
## 
## data:  UsYes and UsNo
## t = 4.9705, df = 354.64, p-value = 1.042e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.7778386 1.7963824
## sample estimates:
## mean of x mean of y 
##  7.866899  6.579789

3.3.1 Interpretación del resultado

Si el p-value del Wilcox test es menor que 0.05 se rechaza la hipótesis inicial y entonces se acepta que existen diferencias significativas entre los dos grupos.

Dado que el número de observaciones es grande, en cualquier caso más de 100 en ambos grupos, los resultados no difieren mucho del t-test.

4 Regresión.

4.1 Modelo de regresión.

Las categorías de la variable ShelveLoc son:

levels(df$ShelveLoc)
## [1] "Bad"    "Good"   "Medium"

Especificamos nivel base para reordenar mediante el método “relevel”, más info: https://www.statology.org/reorder-factor-levels-in-r/

df$ShelveLoc1 = relevel(df$ShelveLoc, ref="Bad")
df$US1 = relevel(df$US, ref="Yes")
df$Urban1 = relevel(df$Urban, ref="Yes")

Evaluamos su correlación para comprobar la colinealidad

round(cor(df_c, method = 'pearson'), 2)  # Una tabla sobre la correlación entre las variables numéricas
##             Sales CompPrice Income Advertising Population Price   Age Education
## Sales        1.00      0.06   0.16        0.31       0.06 -0.43 -0.22     -0.04
## CompPrice    0.06      1.00  -0.08       -0.02      -0.09  0.58 -0.10      0.03
## Income       0.16     -0.08   1.00        0.06      -0.01 -0.06  0.00     -0.06
## Advertising  0.31     -0.02   0.06        1.00       0.27  0.04  0.00     -0.03
## Population   0.06     -0.09  -0.01        0.27       1.00 -0.01 -0.04     -0.11
## Price       -0.43      0.58  -0.06        0.04      -0.01  1.00 -0.10      0.01
## Age         -0.22     -0.10   0.00        0.00      -0.04 -0.10  1.00      0.01
## Education   -0.04      0.03  -0.06       -0.03      -0.11  0.01  0.01      1.00

No existen valores por encima de 0.80. No existen evidencias de multicolinealidad. Comprobamos sus distribuciones.

multi.hist(x = df_c, dcol = c("blue", "red"), dlty = c("dotted", "solid"),
           main = NULL)

Con una baja correlación, existen distribuciones que invitan a la discretización, lo evaluaremos. Mediante una correlación lineal mediante la variable Sales como dependiente.

mod4.1<- lm(Sales~Price+Advertising+Age+Population+ShelveLoc1+US1+Urban1+Education, data=df)
summary(mod4.1)
## 
## Call:
## lm(formula = Sales ~ Price + Advertising + Age + Population + 
##     ShelveLoc1 + US1 + Urban1 + Education, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7366 -1.0133  0.0223  1.0950  4.5818 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      14.1566188  0.7182905  19.709  < 2e-16 ***
## Price            -0.0578371  0.0034677 -16.679  < 2e-16 ***
## Advertising       0.1140733  0.0177141   6.440 3.52e-10 ***
## Age              -0.0469930  0.0050616  -9.284  < 2e-16 ***
## Population       -0.0003609  0.0005879  -0.614    0.540    
## ShelveLoc1Good    4.3905473  0.2436489  18.020  < 2e-16 ***
## ShelveLoc1Medium  1.9421442  0.2004832   9.687  < 2e-16 ***
## US1No            -0.2150244  0.2383507  -0.902    0.368    
## Urban1No         -0.1938727  0.1797412  -1.079    0.281    
## Education        -0.0048699  0.0313835  -0.155    0.877    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.624 on 390 degrees of freedom
## Multiple R-squared:  0.6548, Adjusted R-squared:  0.6469 
## F-statistic: 82.21 on 9 and 390 DF,  p-value: < 2.2e-16

4.2 Interpretar el modelo

Primeras valoraciones:

R-squared->0.6548 significa que puede explicar el 65% de la variabilidad observada en la variable dependiente, “Sales”. Y esto con una significancia ya que el p-value es inferior al 0.05, así que el modelo no es azaroso, aunque un valor inferior a 0.8 nos debe hacer ser cautos. Sobre los regresores:

Hay 4 variables cuya contribución al modelo es poco significativa:

  • Población

  • Us1No -> US

  • Urban1No -> Urban

  • Education

AKAIKE Information Criterion (AIC) es una medida de calidad relativa de un modelo estadístico y proporciona un medio para la selección de un modelo[1]. La función step del paquete stats nos ofrece un algoritmo en el que va probando modelos con el fin de ofrecer lo que mejor puntuación AIC da. En nuestro caso, cuando ejecuto lo que hace es algo previsible, termina eliminando las variables no significativas. Más info: https://www.statology.org/stepwise-regression-r/

slm<-step(object=mod4.1, direction = 'both')
## Start:  AIC=397.7
## Sales ~ Price + Advertising + Age + Population + ShelveLoc1 + 
##     US1 + Urban1 + Education
## 
##               Df Sum of Sq    RSS    AIC
## - Education    1      0.06 1028.4 395.73
## - Population   1      0.99 1029.4 396.09
## - US1          1      2.15 1030.5 396.54
## - Urban1       1      3.07 1031.4 396.90
## <none>                     1028.4 397.70
## - Advertising  1    109.35 1137.7 436.12
## - Age          1    227.28 1255.7 475.58
## - Price        1    733.53 1761.9 611.07
## - ShelveLoc1   2    857.89 1886.3 636.35
## 
## Step:  AIC=395.73
## Sales ~ Price + Advertising + Age + Population + ShelveLoc1 + 
##     US1 + Urban1
## 
##               Df Sum of Sq    RSS    AIC
## - Population   1      0.95 1029.4 394.10
## - US1          1      2.24 1030.7 394.60
## - Urban1       1      3.11 1031.5 394.94
## <none>                     1028.4 395.73
## + Education    1      0.06 1028.4 397.70
## - Advertising  1    109.47 1137.9 434.19
## - Age          1    227.34 1255.8 473.61
## - Price        1    734.01 1762.4 609.19
## - ShelveLoc1   2    859.21 1887.6 634.65
## 
## Step:  AIC=394.1
## Sales ~ Price + Advertising + Age + ShelveLoc1 + US1 + Urban1
## 
##               Df Sum of Sq    RSS    AIC
## - US1          1      2.85 1032.2 393.21
## - Urban1       1      3.35 1032.7 393.40
## <none>                     1029.4 394.10
## + Population   1      0.95 1028.4 395.73
## + Education    1      0.02 1029.4 396.09
## - Advertising  1    114.34 1143.7 434.23
## - Age          1    226.56 1255.9 471.67
## - Price        1    733.29 1762.7 607.25
## - ShelveLoc1   2    862.21 1891.6 633.48
## 
## Step:  AIC=393.21
## Sales ~ Price + Advertising + Age + ShelveLoc1 + Urban1
## 
##               Df Sum of Sq    RSS    AIC
## - Urban1       1      3.48 1035.7 392.55
## <none>                     1032.2 393.21
## + US1          1      2.85 1029.4 394.10
## + Population   1      1.57 1030.7 394.60
## + Education    1      0.07 1032.2 395.18
## - Age          1    225.49 1257.7 470.24
## - Advertising  1    262.99 1295.2 481.99
## - Price        1    730.89 1763.1 605.35
## - ShelveLoc1   2    864.29 1896.5 632.53
## 
## Step:  AIC=392.55
## Sales ~ Price + Advertising + Age + ShelveLoc1
## 
##               Df Sum of Sq    RSS    AIC
## <none>                     1035.7 392.55
## + Urban1       1      3.48 1032.2 393.21
## + US1          1      2.99 1032.7 393.40
## + Population   1      1.89 1033.8 393.82
## + Education    1      0.11 1035.6 394.51
## - Age          1    223.73 1259.4 468.78
## - Advertising  1    266.17 1301.9 482.04
## - Price        1    727.60 1763.3 603.40
## - ShelveLoc1   2    860.98 1896.7 630.56
summary(slm)
## 
## Call:
## lm(formula = Sales ~ Price + Advertising + Age + ShelveLoc1, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6646 -1.0472 -0.0266  1.0462  4.2386 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      13.746428   0.526035  26.132   <2e-16 ***
## Price            -0.057464   0.003454 -16.637   <2e-16 ***
## Advertising       0.123138   0.012237  10.062   <2e-16 ***
## Age              -0.046546   0.005045  -9.225   <2e-16 ***
## ShelveLoc1Good    4.378652   0.242193  18.079   <2e-16 ***
## ShelveLoc1Medium  1.921563   0.198860   9.663   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.621 on 394 degrees of freedom
## Multiple R-squared:  0.6524, Adjusted R-squared:  0.648 
## F-statistic: 147.9 on 5 and 394 DF,  p-value: < 2.2e-16
slm$anova
##           Step Df   Deviance Resid. Df Resid. Dev      AIC
## 1              NA         NA       390   1028.364 397.7042
## 2  - Education  1 0.06349126       391   1028.428 395.7289
## 3 - Population  1 0.94846382       392   1029.376 394.0976
## 4        - US1  1 2.85463700       393   1032.231 393.2053
## 5     - Urban1  1 3.48290617       394   1035.714 392.5527

El modelo no mejora demasiado, puesto que la bondad expresada en R-Squared es muy similar.

4.2.1 modelo con US-no, Urban-no, ShelveLoc- Bad

La expresión matemática del modelo sería:

Salas = 14.09 - 0.0578·Price + 0.1138·Adversitinos - 0.047·age - 0.00035·Population + ShelveLoc=Good(1)·4.3917 + ShelveLoc=Medium(1)·1.94 - US=NO(1)·0.22 - Urban=NO(1)· 0.195 - 0.00486 * Education

Lo que se propone es reconstruir el modelo sabiendo que US=No, Urban=No y ShelveLoc=Bad.

Sustituimos las variables:

OS=No -> 0.22

Urban=No -> 0.195

ShelveLoc=Bad -> ShelveLocGood y ShelveLocMedium = 0

Quedaría:

Salas = 13,676074600 - 0.0578·Price + 0.1138·Adversitinos - 0.047·age - 0.00035·Population

4.3 Predicción

Predecir las ventas con estos datos:

OS=No

Urban = NO

Price = 131

Advertising = 0

población = 139

Age = 40

ShelveLoc = Bad

tenda <- data.frame(Price=131, Advertising=0, Age=40, 139, Population=139, ShelveLoc1="Bad", US1="No", Urban1="No", Education=0)
predict(mod4.1, tenda, type='response')
##        1 
## 4.241176

Sales = 4.169583

Comparamos ahora con una tienda con las características:

tenda1 <- data.frame(Price=131, Advertising=9, Age=40, 139, Population=139, ShelveLoc1="Good", US1="No", Urban1="No", Education=0)
predict(mod4.1, tenda1, type='response')
##        1 
## 9.658382

Sólo ha variado 2 datos del modelo:

Adversiting, que pasa de 0 -> 9.

ShelveLoc que pasa de Bad -> Good.

Según el modelo, por el incremento de la publicidad aumentaría 0.11*9, aprox 1 punto.

De pasar de Bad a Good aumenta los 4.3917379 restantes.

La predicción de Sales pasa de 4.169.583 a 9.586.383 en gran medida por la clase “good”. ## Discretización de la variable Advertising

Probamos el modelo discretizando la variable Advertising. El histograma lo hemos visto antes. Probamos diferentes métodos que nos ofrece el paquete Arules y escogemos la selección de los cortes por clustering en 3 grupos.

table(discretize(df$Advertising, method = "cluster", breaks = 3))
## 
##    [0,2.99) [2.99,9.94)   [9.94,29] 
##         161          89         150
hist(df$Advertising, breaks = 20, main = "K-means")
abline(v = discretize(df$Advertising, method = "cluster", breaks = 3, 
  onlycuts = TRUE), col = "red")

Ejecutamos

df$Advertising_d<-as.factor(discretize(df$Advertising, method = "cluster", breaks = 3, labels = FALSE))

Hemos creado otra variable categórica. Dependiendo de su valor tenemos 3 grupos:

1- [0,5.69) que existen 206 observaciones. 2- [5.69,14.2) con 138. 3- [14.2,29] con 56.

¿Esta discretización nos ofrecerá mejor bondad del modelo de regresión? Comprobémoslo.

mod4.4<- lm(Sales~Price+Advertising_d+Age+Population+ShelveLoc1+US1+Urban1+Education, data=df)
summary(mod4.4)
## 
## Call:
## lm(formula = Sales ~ Price + Advertising_d + Age + Population + 
##     ShelveLoc1 + US1 + Urban1 + Education, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6803 -1.0577 -0.0078  1.0908  4.3216 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.428e+01  7.381e-01  19.343  < 2e-16 ***
## Price            -5.785e-02  3.529e-03 -16.392  < 2e-16 ***
## Advertising_d2    7.567e-01  2.611e-01   2.898  0.00397 ** 
## Advertising_d3    1.553e+00  2.937e-01   5.288 2.07e-07 ***
## Age              -4.792e-02  5.149e-03  -9.306  < 2e-16 ***
## Population       -1.155e-05  5.892e-04  -0.020  0.98437    
## ShelveLoc1Good    4.412e+00  2.476e-01  17.819  < 2e-16 ***
## ShelveLoc1Medium  1.959e+00  2.039e-01   9.608  < 2e-16 ***
## US1No            -4.507e-01  2.564e-01  -1.758  0.07957 .  
## Urban1No         -2.098e-01  1.829e-01  -1.148  0.25186    
## Education         1.844e-03  3.189e-02   0.058  0.95391    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.651 on 389 degrees of freedom
## Multiple R-squared:  0.6442, Adjusted R-squared:  0.635 
## F-statistic: 70.42 on 10 and 389 DF,  p-value: < 2.2e-16

No mejora demasiado el modelo con una R-squared de 0.6432, y vemos que el hecho de la nueva variable Adversiting1 (discretizada) ser del grupo 3 (más gasto) es significativo, menos el nivel 2 de gasto. No discretizamos la variable Advertising*. Como curiosidad, el hecho de hacer un alto gasto en publicidad tendría menos impacto (1.603) que el colocar el poroducto en una estantería “medium” (1.988) y lejos del hecho de utilizarlo “Good” ( 4.399).

5 Análisis de la varianza (ANOVA)

5.1 Anova de un factor

5.1.1 Ventas de calidad de la localización dentro del expositor

ANOVA de significación de la variable ShelveLoc en la variable Sales. #### Visualización

ggplot(df, aes(x=ShelveLoc, y=Sales, color=ShelveLoc))+ geom_boxplot()+theme_bw()

Hay datos atípicos, que ya comentamos en puntos anteriores. Concretamente en la clase “Medium”. La distribución no es simétrica. Parece que las cajas son similares y podría no existir heterocedasticidad.

5.1.1.1 Condiciones aplicación ANOVA. Normalidad.

Los tres subgrupos representan más del 10% (n>40, N>400), siendo independientes entre ellos. Dibujamos la distribución intercuartílica de los grupos.

par(mfrow = c(1,3))
qqnorm(df[df$ShelveLoc == "Bad","Sales"], main = "Bad")
qqline(df[df$ShelveLoc == "Bad","Sales"])
qqnorm(df[df$ShelveLoc == "Good","Sales"], main = "Good")
qqline(df[df$ShelveLoc == "Good","Sales"])
qqnorm(df[df$ShelveLoc == "Medium","Sales"], main = "Medium")
qqline(df[df$ShelveLoc == "Medium","Sales"])

Parece que pasan el test de normalidad, con n>50 utilizamos la corrección de Lilliefors del test de Kolmogorov-Smirnov.

require(nortest) # Este método requiere el paquete nortest
## Loading required package: nortest
by(data = df,INDICES = df$ShelveLoc,FUN = function(x){ lillie.test(x$Sales)})
## df$ShelveLoc: Bad
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x$Sales
## D = 0.079919, p-value = 0.1374
## 
## ------------------------------------------------------------ 
## df$ShelveLoc: Good
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x$Sales
## D = 0.054845, p-value = 0.7634
## 
## ------------------------------------------------------------ 
## df$ShelveLoc: Medium
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x$Sales
## D = 0.037603, p-value = 0.633

Y efectivamente ninguno de los grupos presenta una carencia de normalidad.

5.1.1.2 Condiciones aplicación ANOVA. Homocedasticidad.

Utilizamos el test de Bartlett para evaluar la carencia de homocedasticidad.

bartlett.test(Sales~ShelveLoc, df)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Sales by ShelveLoc
## Bartlett's K-squared = 0.7014, df = 2, p-value = 0.7042

El p-value es muy alto, no existen evidencias de heterocedasticidad.

5.1.1.3 Hipótesis nula y alteranativa

\(H_{0} : \mu_{1}=\mu_{2}\)

\(H_{1} : \mu_{1}\ne\mu_{2}\)

Es decir, como hipótesis nula, las variancias de las poblaciones estudiadas son iguales.

5.1.1.4 Model

anova<- aov(df$Sales ~ df$ShelveLoc)
summary(anova)
##               Df Sum Sq Mean Sq F value Pr(>F)    
## df$ShelveLoc   2  832.8   416.4   77.02 <2e-16 ***
## Residuals    397 2146.5     5.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Con una F de Fisher-Snedecor de 77 nos ofrece una p-value muy inferior al 0.05 podemos afirmar que existe diferencia significiativa entre las variancias de los grupos estudiados.

La suma de los cuadrados entre los grupos es de 832,8.

La suma de los cuadrados total es de 2146.

Las medias cuadráticas vienen dadas por la ratio Sum Sq/Grados libertad(Df) -> 832.8/2 = 416.4 …

F-Value viene dado por la ratio de las medias cuadradas (416.4/5.4=77.02).

Por último el estadístico de Fisher-Snedecor proporciona unas tablas y con un valor de 77 el p-value es inferior a 0.05.

Cómo ha resultado significativa podemos estudiar entre los grupos los resultados, con la función TukeyHSD.

TukeyHSD(anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = df$Sales ~ df$ShelveLoc)
## 
## $`df$ShelveLoc`
##                  diff       lwr       upr p adj
## Good-Bad     4.284730  3.470028  5.099433     0
## Medium-Bad   1.783659  1.114079  2.453239     0
## Medium-Good -2.501072 -3.200125 -1.802019     0

diff -> diferencias entre las medias de ambos grupos. lwr, upr -> intervalo de confianza al 95%. p adj -> p-value después de los ajustes por las comparaciones múltiples. En ese caso, todo significativo.

plot(TukeyHSD(anova))

En el gráfico queda mucho más evidente los intervalos de confianza de los tres grupos, evidenciando la significación de las diferencias entre las medias de las variancias.

5.1.2 Cálculos

  1. Dibujamos gráfico de cajas con orden Bad, Medium, Good.
ggplot(df, aes(x= reorder (ShelveLoc, +Sales), y=Sales, color=ShelveLoc))+ geom_boxplot()+theme_bw()

  1. Cálculo manual ANOVA

Tenemos 3 grupos.

Calculamos tanto las medias y desviación tipo como el número observaciones:

aggregate( Sales ~ ShelveLoc, df, FUN = function(x) c(mean = mean(x), n = length(x), sd = sd(x)))
##   ShelveLoc Sales.mean    Sales.n   Sales.sd
## 1       Bad   5.522917  96.000000   2.356349
## 2      Good   9.807647  85.000000   2.437948
## 3    Medium   7.306575 219.000000   2.266373

La gran media es

mean(df$Sales)
## [1] 7.409975

Msa - Cuadrado de la media entre grups/Grados libertad = 2 (3 grups -1)

Msa = (96*(5.522917 - 7.409975)^2+ 85*(9.807647 - 7.409975)^2 + 219*(7.306575 - 7.409975)^2)/2
Mse = (95*2.356349^2 + 84*2.437948^2 + 218*2.266373^2)/(95+84+218)
F = Msa/Mse
F
## [1] 77.01906
  1. Efectos de los niveles del factor ShelveLoc y la estimación de la varianza del error.

La suma de los cuadrados intra-grupos es el numerador de Msa. La suma de los cuadrados residuales es el numerador de Mse. La suma de los dos anteriores es la suma de cuadrados total.

El valor que se pide es el ratio del 1er/3er.

n2=(96*(5.522917 - 7.409975)^2+ 85*(9.807647 - 7.409975)^2 + 219*(7.306575 - 7.409975)^2)/(Msa*2+Mse*397)
n2
## [1] 0.2795417

Como curiosidad, si queremos calcular este valor, podemos ejecutar la función etaSquared del paquete lsr

require(lsr)
etaSquared(anova)
##                 eta.sq eta.sq.part
## df$ShelveLoc 0.2795417   0.2795417

La varianza del error se calcula mediante el numerador de Mse.

Mse*397
## [1] 2146.483

5.1.3 Interpretación de los resultados

La p-value es inferior a 0.05, por lo que sí existe diferencia significativa entre las medias de los grupos.

Pero, ¿cómo de mayor es esta diferencia? La respuesta nos la puede dar el valor de la estimación de los efectos, en inglés sería eta Squared.

Para conextualizar el valor obtenido antes (0.2795) recurriremos a la literatura de los escritos de Cohen and Miles & Shevin (2001)[2] Diremos que un valor de 0.01 es un valor pequeño, 0.06 medio y más de 0.14 es grande. Ahora si podemos decir que las diferencias son significativas y según la evaluación de la eta squared evaluada por Cohen con 0.2795 es grande.

5.2 Adecuación del modelo.

5.2.1 Normalidad de los residuos.

La normalidad de los residuos la podemos evaluar visualmente mediante el 2º plote del modelo.

plot(anova, 2)

Podemos calcular paramétricamente, aunque visualmente no se aprecia carencia de normalidad.

#para obtener los valores residuales
anova_residuos <- residuals(object = anova)
#Ahora ejecutamos el test de Shapiro-Wilk
shapiro.test(anova_residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  anova_residuos
## W = 0.99574, p-value = 0.3522

Con un p-value superior a 0,05 nos confirma la normalidad de la distribución de los residuos

5.2.2 Homocedasticidad de los residuos

Corresponde al primer plote del modelo.

plot(anova, 1)

En el dibujo no se aprecia heterocedasticidad, pero podemos contrastarlo con el test de levene..

leveneTest(Sales~ShelveLoc, data = df)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.4476 0.6395
##       397

Y se confirma la carencia de heterocedasticidad. Otros gráficos que nos pueden ayudar a interpretar visualmente son:

Scale-Location .- Que nos muestra si los residuales tienen un abanico similar representado por la línea roja que separa los puntos por la mitad de cada grupo.

plot(anova, 3)

El gráfico Residuales vs Leverage nos ayuda a detectar observaciones que puedan influir para que puedan ser detectados como outliers. Corresponde al quinto dibujo del modelo.

plot(anova, 5)

plot(anova, 4)

En este gráfico detectamos las observaciones que pueden influir en función de la distancia de Cook.

Estos puntos ya se detectaron en las primeras box-plot donde nos indicaban que había algunas observaciones fuera de rangos concretamente en el grupo “medium”.

5.3 ANOVA no paramétrico.

Se nos pide evaluar sólo un grupo de observaciones, concretamente aquellas que sean mayores en la mediana de las variables Advertising. Creamos la dataset:

df5_3 <- df[df[,'Advertising']>median(df$Advertising),]

Ahora sólo tenemos 194 observaciones.

Comprobamos criterios ANOVA. Dibujamos los box-plot por grupo:

ggplot(df5_3, aes(x= reorder (ShelveLoc, +Sales), y=Sales, color=ShelveLoc))+ geom_boxplot()+theme_bw()

Estudiamos si pasan test normalidad

require(nortest)
by(data = df5_3, INDICES = df5_3$ShelveLoc, FUN = function(x){ lillie.test(x$Sales)})
## df5_3$ShelveLoc: Bad
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x$Sales
## D = 0.084507, p-value = 0.5987
## 
## ------------------------------------------------------------ 
## df5_3$ShelveLoc: Good
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x$Sales
## D = 0.075913, p-value = 0.728
## 
## ------------------------------------------------------------ 
## df5_3$ShelveLoc: Medium
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x$Sales
## D = 0.077078, p-value = 0.1346

Los tres grupos pasan el test de normalidad.

Utilizamos el test de Bartlett para evaluar la carencia de homocedasticidad.

bartlett.test(Sales~ShelveLoc, df5_3)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Sales by ShelveLoc
## Bartlett's K-squared = 4.2229, df = 2, p-value = 0.1211

No se aprecia heterocedasticidad.

Estudiamos otros parámetros:

anova53<-aov(df5_3$Sales ~ df5_3$ShelveLoc)
summary(anova53)
##                  Df Sum Sq Mean Sq F value Pr(>F)    
## df5_3$ShelveLoc   2  503.9  251.94    49.2 <2e-16 ***
## Residuals       191  978.1    5.12                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hay diferencias significativas. Estudiamos la normalidad y homocedasticidad de las varianzas.

plot(anova53,2)

No se detecta carencia de normalidad.

Evaluamos la heterocedasticidad.

plot(anova53, 1)

anova53_residus <- residuals(object = anova53)
shapiro.test(x = anova53_residus)
## 
##  Shapiro-Wilk normality test
## 
## data:  anova53_residus
## W = 0.99455, p-value = 0.7042

Evaluamos eta Squared

require(lsr)
etaSquared(anova53)
##                    eta.sq eta.sq.part
## df5_3$ShelveLoc 0.3400078   0.3400078

Un valo superior encara que la dataset complerta.

Comparamos 2 a 2.

pairwise.t.test(x = df5_3$Sales, g = df5_3$ShelveLoc, p.adjust.method = "holm", 
    pool.sd = TRUE, paired = FALSE, alternative = "two.sided")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df5_3$Sales and df5_3$ShelveLoc 
## 
##        Bad     Good   
## Good   < 2e-16 -      
## Medium 1.0e-06 7.2e-10
## 
## P value adjustment method: holm

Que confirma la significancia de las medias.

5.3.1 Kruskal test

Ejecutamos el test Kruskal aunque no hemos encontrado violación de los criterios ANOVA. Más info: https://www.statology.org/kruskal-wallis-test-in-r/

kruskal.test(Sales~ShelveLoc, data = df5_3)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Sales by ShelveLoc
## Kruskal-Wallis chi-squared = 64.571, df = 2, p-value = 9.52e-15

Con un p-value inferior a 0.05, podemos concluir que las diferencias entre los grupos son significativas.

La evaluación entre grupos en estadístico no paramétricos sería similar a la función “pairwise.t.test”, pero wilcox.

pairwise.wilcox.test(df5_3$Sales, df5_3$ShelveLoc, p.adjust.method = 'BH')
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df5_3$Sales and df5_3$ShelveLoc 
## 
##        Bad     Good   
## Good   2.1e-11 -      
## Medium 1.4e-05 1.3e-09
## 
## P value adjustment method: BH

Similares datos que confirman la significancia de las diferencias entre los grupos analizados.

6 ANOVA multifactorial

6.1 FActores: ShelveLoc y US

6.1.1 Análisis visual de los efectos principales y posibles interacciones

Gráfico de la variable Salas en función de ShelveLoc y en función de US.

6.1.1.1 Tabla medias variables

Calculamos las medias según los factores por variables

o en taules separades

taula61a<-df %>%
  group_by(ShelveLoc)%>%
  summarize(mean_size = mean(Sales))
taula61b<-df %>%
  group_by(US)%>%
  summarize(mean_size = mean(Sales))

Las mostramos en tablas

print(taula61a)
## # A tibble: 3 × 2
##   ShelveLoc mean_size
##   <fct>         <dbl>
## 1 Bad            5.52
## 2 Good           9.81
## 3 Medium         7.31
print(taula61b)
## # A tibble: 2 × 2
##   US    mean_size
##   <fct>     <dbl>
## 1 No         6.58
## 2 Yes        7.87

6.1.1.2 Gráficos

ggplot(data = df, aes(x = ShelveLoc, y = Sales, colour = US, group = US)) + 
    stat_summary(fun = mean, geom = "point") + stat_summary(fun= mean, 
    geom = "line") + labs(y = "Mitjana (Sales)") + theme_bw()

Otra alternativa sería con la librería ggpubr

ggline(df, x = "ShelveLoc", y = "Sales", color = "US",
       add = c("mean", "jitter"),
       palette = c("#00AFBB", "#E7B800"))

La media de ventas, independientemente de la situación de la variable ShelveLoc es superior, y mayor cuando la calidad aumenta. Los puntos de colores se encuentran muy mezclados en “Medium” que es donde más observaciones existen. En cambio en “good” parece que hay mayor separación de los colores, discriminaría mejor.

7 Modelo ANOVA two-way

Más información sobre anovas https://cienciadedatos.net/documentos/19_anova

anova2 <- aov(Sales~ShelveLoc * US, data= df)
summary(anova2)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## ShelveLoc      2  832.8   416.4  83.764  < 2e-16 ***
## US             1  115.7   115.7  23.272 2.01e-06 ***
## ShelveLoc:US   2   72.1    36.0   7.248 0.000811 ***
## Residuals    394 1958.7     5.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nos confirma lo que visualmente se apreciaba, con valores significativos. Hay correlación por separado, tanto ShelveLoc como US tienen una influencia en Salas individualmente, pero también de forma conjunta entre ellos. Evaluamos el efecto

etaSquared(anova2)
##                  eta.sq eta.sq.part
## ShelveLoc    0.26744486  0.28916585
## US           0.03883167  0.05577093
## ShelveLoc:US 0.02418750  0.03548497

Aquí ya podemos ver por separado que la variable ShelveLoc aporta un efecto grande, siendo medio tanto OS como la combinación de los dos. Tal y como hemos hecho antes, evaluamos los gráficos de los residuos:

plot(anova2, 1)

Indica que existen 3 puntos que pueden influir, son los ya señalados y que se sugeriría eliminar. Pasamos el test de homocedasticidad de los residuos con el leveneTest

leveneTest(Sales~ShelveLoc*US, data = df)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   5  1.4632 0.2009
##       394

No hay indicios de carencia de homocedasticidad con una p >0.05

plot(anova2, 2)

Indica las 3 observaciones anteriores. confirmamos con shapiro

anova2_residus <- residuals(object=anova2)
shapiro.test(anova2_residus)
## 
##  Shapiro-Wilk normality test
## 
## data:  anova2_residus
## W = 0.99616, p-value = 0.4463

No existen indicios de falta de normalidad en la variancia de los residuos.

He detectado una carencia de uniformidad en las distribuciones de las clases. 1 Bad No 34 5.29 2 Bad Yes 62 5.65 3 Good No 24 7.71 4 Good Yes 61 10.6 5 Medium No 84 6.78 6 Medium Yes 135 7.63

En este caso se puede evaluar para ver si concide el resultado mediante la función Anova del paquete car[3]. El tipo III se puede evaluar por cuándo hay efecto e interacción

Anova(anova2, type='III')
## Anova Table (Type III tests)
## 
## Response: Sales
##               Sum Sq  Df  F value    Pr(>F)    
## (Intercept)   950.30   1 191.1529 < 2.2e-16 ***
## ShelveLoc      90.91   2   9.1438 0.0001313 ***
## US              2.94   1   0.5906 0.4426632    
## ShelveLoc:US   72.06   2   7.2477 0.0008107 ***
## Residuals    1958.73 394                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Este test, que tendría en cuenta la falta de observaciones similares entre los grupos, nos dice que OS ya no es significiativo. Este resultado nos indica que deberemos ser conservadores a la hora de dar por sentado que US influye de manera significativa en las ventas, dado la falta de similitud entre el número de observaciones.

Podemos trabajar con mayor profundidad esta relación con otros esos. Hagamos lo siguiente:

lmAnova = lm(Sales~ShelveLoc + US, data=df)

Para trabajar el siguiente gráfico utilizo el paquete phia.

IM <- interactionMeans(lmAnova)
plot(IM)

En estos gráficos podría interpretarse como que la variabilidad correspondería en mayor medida a la variable ShelveLoc que a US.

7.1 Factores ShelveLoc y Urban

7.1.1 Análisis visual de los efectos principales y posibles interacciones.

Haremos una sola mesa

taula62<-group_by(df, ShelveLoc, Urban)%>%
  summarise(
    count= n(),
    mean= mean(Sales)
  )
## `summarise()` has grouped output by 'ShelveLoc'. You can override using the
## `.groups` argument.
print(taula62)
## # A tibble: 6 × 4
## # Groups:   ShelveLoc [3]
##   ShelveLoc Urban count  mean
##   <fct>     <fct> <int> <dbl>
## 1 Bad       No       22  5.55
## 2 Bad       Yes      74  5.52
## 3 Good      No       28  9.70
## 4 Good      Yes      57  9.86
## 5 Medium    No       68  7.24
## 6 Medium    Yes     151  7.34

Dibujaremos los gráficos

ggline(df, x = "ShelveLoc", y = "Sales", color = "Urban",
       add = c("mean", "jitter"),
       palette = c("#00AFBB", "#E7B800"))

Se aprecia gran mezcla de colores por lo que costaría discriminar por Urban, veremos los datos que nos indican.

Ejecutamos el modelo

anova3 <- aov(Sales~ShelveLoc * Urban, data= df)
summary(anova3)
##                  Df Sum Sq Mean Sq F value Pr(>F)    
## ShelveLoc         2  832.8   416.4  76.469 <2e-16 ***
## Urban             1    0.6     0.6   0.105  0.746    
## ShelveLoc:Urban   2    0.3     0.2   0.030  0.970    
## Residuals       394 2145.6     5.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Efectivamente, sólo ShelveLoc muestra una diferenciación. El resto con valores muy altos tanto en Urban como la interacción entre ambas variables no indica significación.

etaSquared(anova3)
##                       eta.sq  eta.sq.part
## ShelveLoc       0.2791802253 0.2793659969
## Urban           0.0001920391 0.0002665926
## ShelveLoc:Urban 0.0001114437 0.0001547257

Es un valor similar al que aportaba con ANOVA one-way, Urban no aporta nada al modelo.

leveneTest(Sales~ShelveLoc*Urban, data = df)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   5  0.9589 0.4428
##       394
plot(anova3, 1)

No se aprecia signos de heterogeneidad.

plot(anova3, 2)

anova3_residus <- residuals(object=anova3)
shapiro.test(anova3_residus)
## 
##  Shapiro-Wilk normality test
## 
## data:  anova3_residus
## W = 0.99576, p-value = 0.3563

No se aprecia carencia en la normalidad de la distribución de los residuos. Al igual que antes, no existe una proporcionalidad en las observaciones por grupo, ejecutamos el test Anova, en este caso, tipo 2.

Anova(anova3, type='II')
## Anova Table (Type II tests)
## 
## Response: Sales
##                  Sum Sq  Df F value Pr(>F)    
## ShelveLoc        831.77   2 76.3704 <2e-16 ***
## Urban              0.57   1  0.1051  0.746    
## ShelveLoc:Urban    0.33   2  0.0305  0.970    
## Residuals       2145.58 394                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Confirma totalmente el resultado.

7.1.2 Conclusiones

La comparativa de dos vías con los factores ShelveLoc y Urban: - Pasa los test de normalidad y homogeneidad de los residuos. - El resultado indica que el factor Urban no aporta diferenciación en las medias. - La interacción de ShelveLoc y Urban tampoco demuestran siginificancia en la diferencia de sus medias.

8 Comparaciones múltiples

Utilizamos el scheffe test del peste agricolae para las 2 combinaciones[4]. Si ejecutamos el test con group=TRUE, ofrece esta información.

comp_scheffe<-scheffe.test(anova2, "ShelveLoc", group=TRUE)
comp_scheffe
## $statistics
##    MSerror  Df        F     Mean    CV
##   4.971391 394 3.018626 7.409975 30.09
## 
## $parameters
##      test    name.t ntr alpha
##   Scheffe ShelveLoc   3  0.05
## 
## $means
##           Sales      std   r        se  Min   Max    Q25  Q50     Q75
## Bad    5.522917 2.356349  96 0.2275639 0.37 11.67 4.0525 5.21  7.4625
## Good   9.807647 2.437948  85 0.2418408 3.58 16.27 7.8000 9.54 11.6200
## Medium 7.306575 2.266373 219 0.1506666 0.00 13.36 5.6250 7.38  8.7750
## 
## $comparison
## NULL
## 
## $groups
##           Sales groups
## Good   9.807647      a
## Medium 7.306575      b
## Bad    5.522917      c
## 
## attr(,"class")
## [1] "group"

Pero realmente lo que nos interesa es si hay significación entre los grupos, y por eso evaluamos si las letras en los grupos son distintas:

comp_scheffe$groups
##           Sales groups
## Good   9.807647      a
## Medium 7.306575      b
## Bad    5.522917      c

Y aquí se ve cómo las letras son diferentes, entonces la diferencia entre los grupos es significativa. Si utilizamos group=FALSE en la parte $groups aparecerá NULL en cambio nos mostrará el p-value en $comparison. Como indica, es significativo.

comp_scheffe2<-scheffe.test(anova2, "US",  group = FALSE, console = TRUE)
## 
## Study: anova2 ~ "US"
## 
## Scheffe Test for Sales 
## 
## Mean Square Error  : 4.971391 
## 
## US,  means
## 
##        Sales      std   r        se  Min   Max    Q25  Q50    Q75
## No  6.579789 2.228414 142 0.1871091 0.00 12.01 5.0800 6.61 8.1650
## Yes 7.866899 2.877131 258 0.1388127 0.37 16.27 5.7625 7.79 9.9875
## 
## Alpha: 0.05 ; DF Error: 394 
## Critical Value of F: 3.865169 
## 
## Comparison between treatments means
## 
##          Difference pvalue sig       LCL       UCL
## No - Yes   -1.28711      0 *** -1.745146 -0.829075
comp_scheffe2
## $statistics
##    MSerror  Df        F     Mean    CV
##   4.971391 394 3.865169 7.409975 30.09
## 
## $parameters
##      test name.t ntr alpha
##   Scheffe     US   2  0.05
## 
## $means
##        Sales      std   r        se  Min   Max    Q25  Q50    Q75
## No  6.579789 2.228414 142 0.1871091 0.00 12.01 5.0800 6.61 8.1650
## Yes 7.866899 2.877131 258 0.1388127 0.37 16.27 5.7625 7.79 9.9875
## 
## $comparison
##          Difference pvalue sig       LCL       UCL
## No - Yes   -1.28711      0 *** -1.745146 -0.829075
## 
## $groups
## NULL
## 
## attr(,"class")
## [1] "group"

También podemos utilizar este test post-hoc donde indica que:

Bad:Yes-Bad:No

Medium:No-Good:No

Medium:Yes-Good:No

Medium:Yes-Medium:No

Estas combinaciones No son significativas. El resto sí.

TukeyHSD(anova2)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sales ~ ShelveLoc * US, data = df)
## 
## $ShelveLoc
##                  diff       lwr       upr p adj
## Good-Bad     4.284730  3.503495  5.065966     0
## Medium-Bad   1.783659  1.141584  2.425733     0
## Medium-Good -2.501072 -3.171409 -1.830735     0
## 
## $US
##            diff       lwr      upr   p adj
## Yes-No 1.120046 0.6620106 1.578082 2.2e-06
## 
## $`ShelveLoc:US`
##                             diff         lwr         upr     p adj
## Good:No-Bad:No        2.42406863  0.72171089  4.12642637 0.0007766
## Medium:No-Bad:No      1.49323529  0.19532758  2.79114301 0.0136206
## Bad:Yes-Bad:No        0.36565465 -0.99698878  1.72829808 0.9726282
## Good:Yes-Bad:No       5.34585824  3.97926478  6.71245171 0.0000000
## Medium:Yes-Bad:No     2.34745752  1.12222385  3.57269118 0.0000011
## Medium:No-Good:No    -0.93083333 -2.40874517  0.54707850 0.4645112
## Bad:Yes-Good:No      -2.05841398 -3.59348924 -0.52333872 0.0019764
## Good:Yes-Good:No      2.92178962  1.38320694  4.46037229 0.0000014
## Medium:Yes-Good:No   -0.07661111 -1.49112746  1.33790523 0.9999875
## Bad:Yes-Medium:No    -1.12758065 -2.19669178 -0.05846951 0.0319343
## Good:Yes-Medium:No    3.85262295  2.77848180  4.92676410 0.0000000
## Medium:Yes-Medium:No  0.85422222 -0.03313286  1.74157731 0.0667951
## Good:Yes-Bad:Yes      4.98020360  3.82867771  6.13172948 0.0000000
## Medium:Yes-Bad:Yes    1.98180287  1.00219493  2.96141081 0.0000002
## Medium:Yes-Good:Yes  -2.99840073 -3.98349580 -2.01330566 0.0000000

Graficamos

plot(TukeyHSD(anova2))

Otra prueba de que se utiliza, a tenor del encontrado en la literatura y tutoriales en internet es la Prueba de Bonferroni.

sld_comp<-LSD.test(anova2, "ShelveLoc", p.adj='bon', console = TRUE)
## 
## Study: anova2 ~ "ShelveLoc"
## 
## LSD t Test for Sales 
## P value adjustment method: bonferroni 
## 
## Mean Square Error:  4.971391 
## 
## ShelveLoc,  means and individual ( 95 %) CI
## 
##           Sales      std   r        se      LCL       UCL  Min   Max    Q25
## Bad    5.522917 2.356349  96 0.2275639 5.075525  5.970308 0.37 11.67 4.0525
## Good   9.807647 2.437948  85 0.2418408 9.332187 10.283107 3.58 16.27 7.8000
## Medium 7.306575 2.266373 219 0.1506666 7.010364  7.602786 0.00 13.36 5.6250
##         Q50     Q75
## Bad    5.21  7.4625
## Good   9.54 11.6200
## Medium 7.38  8.7750
## 
## Alpha: 0.05 ; DF Error: 394
## Critical Value of t: 2.404246 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##           Sales groups
## Good   9.807647      a
## Medium 7.306575      b
## Bad    5.522917      c
sld_comp
## $statistics
##    MSerror  Df     Mean    CV
##   4.971391 394 7.409975 30.09
## 
## $parameters
##         test  p.ajusted    name.t ntr alpha
##   Fisher-LSD bonferroni ShelveLoc   3  0.05
## 
## $means
##           Sales      std   r        se      LCL       UCL  Min   Max    Q25
## Bad    5.522917 2.356349  96 0.2275639 5.075525  5.970308 0.37 11.67 4.0525
## Good   9.807647 2.437948  85 0.2418408 9.332187 10.283107 3.58 16.27 7.8000
## Medium 7.306575 2.266373 219 0.1506666 7.010364  7.602786 0.00 13.36 5.6250
##         Q50     Q75
## Bad    5.21  7.4625
## Good   9.54 11.6200
## Medium 7.38  8.7750
## 
## $comparison
## NULL
## 
## $groups
##           Sales groups
## Good   9.807647      a
## Medium 7.306575      b
## Bad    5.522917      c
## 
## attr(,"class")
## [1] "group"

Como indica “Treatments with the same letter are not significantly different.” entonces cada grupo tiene una letra distinta.

sld_comp2<-LSD.test(anova2, "US", p.adj='bon', console = TRUE)
## 
## Study: anova2 ~ "US"
## 
## LSD t Test for Sales 
## P value adjustment method: bonferroni 
## 
## Mean Square Error:  4.971391 
## 
## US,  means and individual ( 95 %) CI
## 
##        Sales      std   r        se      LCL      UCL  Min   Max    Q25  Q50
## No  6.579789 2.228414 142 0.1871091 6.211932 6.947646 0.00 12.01 5.0800 6.61
## Yes 7.866899 2.877131 258 0.1388127 7.593993 8.139805 0.37 16.27 5.7625 7.79
##        Q75
## No  8.1650
## Yes 9.9875
## 
## Alpha: 0.05 ; DF Error: 394
## Critical Value of t: 1.966003 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##        Sales groups
## Yes 7.866899      a
## No  6.579789      b
sld_comp2
## $statistics
##    MSerror  Df     Mean    CV
##   4.971391 394 7.409975 30.09
## 
## $parameters
##         test  p.ajusted name.t ntr alpha
##   Fisher-LSD bonferroni     US   2  0.05
## 
## $means
##        Sales      std   r        se      LCL      UCL  Min   Max    Q25  Q50
## No  6.579789 2.228414 142 0.1871091 6.211932 6.947646 0.00 12.01 5.0800 6.61
## Yes 7.866899 2.877131 258 0.1388127 7.593993 8.139805 0.37 16.27 5.7625 7.79
##        Q75
## No  8.1650
## Yes 9.9875
## 
## $comparison
## NULL
## 
## $groups
##        Sales groups
## Yes 7.866899      a
## No  6.579789      b
## 
## attr(,"class")
## [1] "group"

Y también sale significativo por el factor US.

Para comparar, ejecutaré el test de scheffe por el factor URBAN

comp_scheffe3<-scheffe.test(anova3, "Urban",  group = TRUE, console = TRUE)
## 
## Study: anova3 ~ "Urban"
## 
## Scheffe Test for Sales 
## 
## Mean Square Error  : 5.445632 
## 
## Urban,  means
## 
##        Sales      std   r        se  Min   Max   Q25  Q50   Q75
## No  7.509237 2.723436 118 0.2148242 0.00 13.91 5.440 7.67 9.320
## Yes 7.368440 2.740159 282 0.1389631 0.37 16.27 5.375 7.40 9.125
## 
## Alpha: 0.05 ; DF Error: 394 
## Critical Value of F: 3.865169 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Means with the same letter are not significantly different.
## 
##        Sales groups
## No  7.509237      a
## Yes 7.368440      a
comp_scheffe3
## $statistics
##    MSerror  Df        F     Mean       CV
##   5.445632 394 3.865169 7.409975 31.49252
## 
## $parameters
##      test name.t ntr alpha
##   Scheffe  Urban   2  0.05
## 
## $means
##        Sales      std   r        se  Min   Max   Q25  Q50   Q75
## No  7.509237 2.723436 118 0.2148242 0.00 13.91 5.440 7.67 9.320
## Yes 7.368440 2.740159 282 0.1389631 0.37 16.27 5.375 7.40 9.125
## 
## $comparison
## NULL
## 
## $groups
##        Sales groups
## No  7.509237      a
## Yes 7.368440      a
## 
## attr(,"class")
## [1] "group"

Y comprobamos cómo los grupos (Yes y No) tienen la misma letra, y no son significativamente diferentes.

9 Conclusiones

9.1 Sobre la base de datos

Hemos trabajado una base de datos sobre 400 observaciones que provienen de una simulación con 14 variables, de las cuales 3 son categóricas. La variable dependiente a estudiar ha sido las ventas “Salas” y las variables factoriales estudiadas han sido “US”, “Urban” y “ShelveLoc”.

9.2 Análisis gráfico

He presentado un estudio gráfico inicial para visualizar las distribuciones de las variables y he detectado 3 outliers que en otras pruebas se han reiterado, pero no se han tratado. He tomado la decisión de no eliminar las observaciones.

9.3 Contraste medias

He presentado el contraste de medias de la variable US y ha salido significativamente diferente, y esto se ha repetido en los estudios factoriales posteriores. También he demostrado que la variable OS presenta heterocedasticidad y debe compararse mediante test no paramétricos.

9.4 Regresión

En la regresión lineal con la variable dependiendo “Salas” la variable ShelveLoc se presenta como muy influyente con un coeficiente muy importante en la clase “good”. Se presentan varios modelos pero no mejoramos la bondad expresada en \(R^2\). He buscado formas para mejorar el modelo mediante la discretización de la variable Advertising, pero no ha dado el resultado esperado.

9.5 Estudio factorial

Una vez realizado tanto una anova de una vía como de dos vías, como variable dependiendo “Salas” y las variables factoriales, “ShelveLoc”, “Urban” y “US” se aprecia una diferencia significativa entre las medias de ShelveLoc y US . La variable Urban no presenta diferencia significativa.

9.6 Relación resultados ANOVA y Regresión lineal

De las variables estudiadas, las más influyentes en la regresión han sido: Price - una relación negativa respecto al precio de la competencia, lógico pero con poca variabilidad. Advertising - Es la variable continua más inluyente en la regresión lineal.

De las variables factoriales ShelveLoc es la más influyente con mucha diferencia. Si comparamos con la variable eta Squared del modelo anova, se confirma este hecho. Una conclusión es que las variables con mayor significancia en un modelo anova podrían ser las más influyentes en un modelo de regresión lineal, respecto a los otros factores, y además podrían ser las más discriminantes.