Datos3 <- data.frame("Ventas"=c(12,10,15,19,11,11,17,16,14,15,27,
33,22,26,28,23,20,18,17),
"Empaque"=as.factor(c("1","1","1","1","1",
"2","2","2","2","2",
"3","3","3","3","3",
"4","4","4","4")))
boxplot(Ventas ~ Empaque, data=Datos3, col="yellow", outline=FALSE)
De esta gráfica podemos ver que en promedio las ventas del empaque 3 son mayores con respecto a las otras 3 (Empaque 1, Empaque 2, Empaque 4)
levels(Datos3$Empaque)
## [1] "1" "2" "3" "4"
Cambiemos la variable de referencia de la variable tipo factor(empaque) Hagamos el cambio tal que nuestra referencia sea el empaque 4
Datos3$Empaque <- factor(Datos3$Empaque,
levels=levels(Datos3$Empaque)[c(4,1,2,3)])
levels(Datos3$Empaque)
## [1] "4" "1" "2" "3"
Como la variable empaque tiene 4 niveles, a nuestro modelo solo van a entrar 3, la que hace referencia al empaque 1, empaque 2 y empaque 3:
\[ E(Y;X) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_3 \] Donde Y = ventas, X = empaque
\[\begin{align*} &E(Y;X1)=\beta_0 + \beta_1\\ &E(Y;X2)=\beta_0 + \beta_2\\ &E(Y;X3)=\beta_0 + \beta_3\\ &E(Y;X4)=\beta_0 \end{align*}\]
Vamos a ajustar el modelo
# b0 = 19.500
# b1 = -6.100
# b2 = -4.900
# b3 = 7.700
fit <- lm(Ventas ~ Empaque, data = Datos3)
summary(fit)
##
## Call:
## lm(formula = Ventas ~ Empaque, data = Datos3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.20 -1.95 -0.20 1.50 5.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.50 1.62 12.01 4.3e-09 ***
## Empaque1 -6.10 2.18 -2.80 0.013 *
## Empaque2 -4.90 2.18 -2.25 0.040 *
## Empaque3 7.70 2.18 3.53 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.2 on 15 degrees of freedom
## Multiple R-squared: 0.788, Adjusted R-squared: 0.746
## F-statistic: 18.6 on 3 and 15 DF, p-value: 2.58e-05
lo anterior nos queda estimado como
\[\begin{align*} &E(Y;X1)=\beta_0 + \beta_1=19.500 - 6.100 = 13.4\\ &E(Y;X2)=\beta_0 + \beta_2=19.500 - 4.900 = 14.6\\ &E(Y;X3)=\beta_0 + \beta_3= 19.500 + 7.700 = 27.2\\ &E(Y;X4)=\beta_0 = 19.500 \end{align*}\]
En nuestro problema las hipotesis que se contrastan son
\[ H_0: \beta_1=0 \ \ y\ \ \beta_2=0 \ \ y\ \ \beta_3=0 \ \ vs \ \ H_a: \beta_1\neq0\ \ o\ \ \beta_2\neq0 \ \ o\ \ \beta_3\neq0 \]
summary(fit)
##
## Call:
## lm(formula = Ventas ~ Empaque, data = Datos3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.20 -1.95 -0.20 1.50 5.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.50 1.62 12.01 4.3e-09 ***
## Empaque1 -6.10 2.18 -2.80 0.013 *
## Empaque2 -4.90 2.18 -2.25 0.040 *
## Empaque3 7.70 2.18 3.53 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.2 on 15 degrees of freedom
## Multiple R-squared: 0.788, Adjusted R-squared: 0.746
## F-statistic: 18.6 on 3 and 15 DF, p-value: 2.58e-05
Con un \(p-value\) de 2.585e-05 se rechaza nuestra hipotesis \(H_o\). Lo que esto quiere decir es que al menos un empaque tiene un promedio distinto respecto a las ventas contra el nivel de referencia, en este caso contra el diseño 4.
Para contestar esta pregunta vamos a contrastar las siguientes hipótesis
\[ H_0: \beta_1=0\ \ y\ \ \beta_2=0\ \ y\ \ \beta_3=0\ \ vs\ \ H_a: \beta_1\neq0 \ \ o\ \ \beta_2\neq0 \ \ o\ \ \beta_3\neq0 \]
# Haciendo una prueba lineal general
library(multcomp)
K <- matrix(c(0,1,0,0,
0,0,1,0,
0,0,0,1), ncol=4, nrow=3, byrow=TRUE)
m <- c(0,0,0)
summary(glht(fit, linfct=K, rhs=m), test=Ftest())
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 1 == 0 -6.1
## 2 == 0 -4.9
## 3 == 0 7.7
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 18.6 3 15 2.58e-05
Pruebas de hipótesis simultáneas (haciéndolas de par en par) usando tukey:
## Usando Tukey obtenemos lo siguiente
summary(glht(fit, linfct = mcp(Empaque = "Tukey")))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = Ventas ~ Empaque, data = Datos3)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 - 4 == 0 -6.10 2.18 -2.80 0.058 .
## 2 - 4 == 0 -4.90 2.18 -2.25 0.155
## 3 - 4 == 0 7.70 2.18 3.53 0.014 *
## 2 - 1 == 0 1.20 2.05 0.58 0.935
## 3 - 1 == 0 13.80 2.05 6.72 <0.001 ***
## 3 - 2 == 0 12.60 2.05 6.13 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Tomando \(\alpha=0.5\), para cada par tenemos lo siguiente:
De aquí notemos que para cualquier prueba simultánea que involucre a 3, siempre nos dice que “se puede considerar una diferencia de forma simultanea”. Lo cual tiene sentido ya que nuestro empaque 3 es el que se veía que tenía en promedio una esperanza distinta con respecto a los otros empaques
\(H_0\):
\[\begin{align*} &E(Y;X3)=\beta_0 + \beta_3 \leq E(Y;X1)= \beta_0 + \beta_1 \ \ o\\ &E(Y;X3)=\beta_0 + \beta_3 \leq E(Y;X2)= \beta_0 + \beta_2 \ \ o\\ &E(Y;X3)=\beta_0 + \beta_3 \leq E(Y;X4)= \beta_0 \end{align*}\]
versus
\(H_a\):
\[\begin{align*} &E(Y;X3)=\beta_0 + \beta_3 > E(Y;X1)= \beta_0 + \beta_1 \ \ y\\ &E(Y;X3)=\beta_0 + \beta_3 > E(Y;X2)= \beta_0 + \beta_2 \ \ y\\ &E(Y;X3)=\beta_0 + \beta_3 > E(Y;X4)= \beta_0 \end{align*}\]
# Haciendo la prueba
K1 <- matrix(c(0,-1,0,1,
0,0,-1,1,
0,0,0,1), ncol=4, nrow=3, byrow=TRUE)
m1 <- c(0,0,0)
summary(glht(fit, linfct=K1, rhs=m1, alternative="greater"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = Ventas ~ Empaque, data = Datos3)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>t)
## 1 <= 0 13.80 2.05 6.72 <0.001 ***
## 2 <= 0 12.60 2.05 6.13 <0.001 ***
## 3 <= 0 7.70 2.18 3.53 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Notemos que con p-values igual a 0.001 , 0.001, 0.00387 se rechaza la hipótesis nula, en términos de nuestro problema la interpretación es que podemos decir es que el diseño del empaque 3 es el que aumenta las ventas en comparación con el resto de empaques.