knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

library(tidyverse)
## ── Attaching packages ────────────────────
## ✔ ggplot2 3.0.0     ✔ purrr   0.2.5
## ✔ tibble  1.4.2     ✔ dplyr   0.7.6
## ✔ tidyr   0.8.1     ✔ stringr 1.3.1
## ✔ readr   1.1.1     ✔ forcats 0.3.0
## ── Conflicts ──── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library(readxl)
library(kableExtra)
library(jmv)
## 
## Attaching package: 'jmv'
## The following object is masked from 'package:stats':
## 
##     anova
datos <- read_csv(file = "kenton.csv") %>%
  transmute(Local = as.factor(Local),
            Diseño = as.factor(Diseno),
            Venta)

# mean(datos$Venta)

Estos datos corresponden a las ventas de un producto (variable Ventas) según el diseño del empaque (variable Diseno). Por lo tanto, Diseno es un factor con 4 niveles (1, 2, 3, 4), y para cada nivel hay 5 réplicas (variable Local), excepto para el diseño 3 que tiene sólo 4 réplicas.

Podemos ver la cantidad de observaciones por tratamiento mediante

Resumen de observaciones

datos %>% 
  group_by(Diseño) %>%
  summarise(N = n()) %>% 
  kable() %>% 
  kable_styling(
    bootstrap_options = c(
      "striped", 
      "hover"),
    full_width = FALSE
  )
Diseño N
1 5
2 5
3 4
4 5 
# table(datos$Diseño)

Un cajón con bigotes nos puede dar una primera impresión de la distribución de la respuesta para cada nivel del factor.

datos %>% 
  ggplot(aes(x = Diseño, y = Venta)) +
  geom_boxplot()

Algo de estadística descriptiva:

datos %>% 
  group_by(Diseño) %>% 
  summarise(N = n(),
            Promedio = mean(Venta, na.rm = TRUE),
            D.E. = sd(Venta, na.rm = TRUE),
            Min = min(Venta, na.rm = TRUE),
            Max = max(Venta, na.rm = TRUE)) %>% 
  kable() %>% 
  kable_styling(
    bootstrap_options = c(
      "striped", 
      "hover"),
    full_width = FALSE
  )
Diseño N Promedio D.E. Min Max
1 5 14.6 2.302173 11 17
2 5 13.4 3.646917 10 19
3 4 19.5 2.645751 17 23
4 5 27.2 3.962323 22 33

Podemos mejorar un poco la tabla

redondear <- function(x, dígitos = 2) 
{
  format(round(x, digits = dígitos), nsmall = dígitos)
}

media <- function(x, dígitos = 2, na.rm = TRUE)
{
  redondear(mean(x, na.rm = na.rm), dígitos = dígitos)
}

d.e. <- function(x, dígitos = 2, na.rm = TRUE)
{
  redondear(sd(x, na.rm = na.rm), dígitos = dígitos)
}

mínimo <- function(x, dígitos = 2, na.rm = TRUE)
{
  redondear(min(x, na.rm = na.rm), dígitos = dígitos)
}

máximo <- function(x, dígitos = 2, na.rm = TRUE)
{
  redondear(max(x, na.rm = na.rm), dígitos = dígitos)
}

resumen <- datos %>% 
  group_by(Diseño) %>% 
  summarise(N = n(),
            Promedio = media(Venta, na.rm = TRUE, dígitos = 1),
            D.E. = d.e.(Venta, na.rm = TRUE, dígitos = 1),
            Min = mínimo(Venta, na.rm = TRUE, dígitos = 1),
            Max = máximo(Venta, na.rm = TRUE, dígitos = 1)) 

resumen%>% 
  kable() %>% 
  kable_styling(
    bootstrap_options = c(
      "striped", 
      "hover"),
    full_width = FALSE
  )
Diseño N Promedio D.E. Min Max
1 5 14.6 2.3 11.0 17.0
2 5 13.4 3.6 10.0 19.0
3 4 19.5 2.6 17.0 23.0
4 5 27.2 4.0 22.0 33.0

Podemos ver un gráfico más elaborado con la intención de ajustar un modelo ANOVA

datos %>% 
  ggplot(aes(x = Diseño, 
             y = Venta)) +
  geom_point() +
  geom_point(data = resumen,
             aes(x = as.numeric(Diseño),
                 y = as.numeric(Promedio)),
             shape = 4,
             size = 5) +
  geom_line(data = resumen,
             aes(x = as.numeric(Diseño),
                 y = as.numeric(Promedio)),
            size = 1,
            linetype = "dashed",
            alpha = 0.5) +
  geom_hline(yintercept = mean(datos$Venta),
             linetype = "dashed",
             alpha = 0.5)

resumen %>% 
  ggplot(aes(x = as.numeric(Diseño),
             y = as.numeric(Promedio))) +
  geom_point(shape = 4,
             size = 5) +
  geom_line(size = 1,
            linetype = "dashed",
            alpha = 0.5) +
  geom_errorbar(aes(ymin = (as.numeric(Promedio) - as.numeric(D.E.)),
                     ymax = (as.numeric(Promedio) + as.numeric(D.E.))),
                width = 0.1) +
  geom_jitter(data = datos,
              mapping = aes(x = as.numeric(Diseño),
                            y = Venta),
              width = 0.1) +
  xlab("Diseño") +
  ylab("Ventas")

df_shapes <- data.frame(shape = 0:24)
ggplot(df_shapes, aes(0, 0, shape = shape)) +
  geom_point(aes(shape = shape), size = 5, fill = 'red') +
  scale_shape_identity() +
  facet_wrap(~shape) +
  theme_void()

Ajuste de modelo ANOVA

modelo1 <- lm(Venta ~ Diseño, data = datos)

# modelo1 <- datos %>% lm(Venta ~ Diseño, data = .)

Resumen

summary(modelo1)
## 
## Call:
## lm(formula = Venta ~ Diseño, data = datos)
## 
## 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)   14.600      1.452  10.053 4.66e-08 ***
## Diseño2       -1.200      2.054  -0.584   0.5677    
## Diseño3        4.900      2.179   2.249   0.0399 *  
## Diseño4       12.600      2.054   6.135 1.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.248 on 15 degrees of freedom
## Multiple R-squared:  0.7881, Adjusted R-squared:  0.7457 
## F-statistic: 18.59 on 3 and 15 DF,  p-value: 2.585e-05

Tabla ANOVA

stats::anova(modelo1)
## Analysis of Variance Table
## 
## Response: Venta
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Diseño     3 588.22 196.074  18.591 2.585e-05 ***
## Residuals 15 158.20  10.547                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
588.22/3
## [1] 196.0733
158.2/15
## [1] 10.54667
196.07/10.547
## [1] 18.59012
tibble(x = seq(0, 25, 0.5),
       y = df(x = x, df1 = 3, df2 = 15)) %>% 
  ggplot(aes(x = x, y = y)) + 
  geom_line() +
  geom_vline(xintercept = 18.591,
             linetype = "dashed",
             alpha = 0.5)

1 - pf(q = 18.591, df1 = 3, df2 = 15)
## [1] 2.585007e-05

Análisis de la varianza

La variación total corresponde a las desviaciones de las observaciones $ Y_i $ con respecto a la media general $ {Y} $. La suma de estas desviaciones al cuadrado se mide con la Suma de Cuadrados Totales \[SSTO=\sum{(Y_i-\bar{Y})^2}\]

Cuando tomamos en consideración las covariables $ X $, podemos medir cómo se aproxima el ajuste lineal $ $ al valor observado $ Y $. Las desviaciones entre el valor ajustado y el observado son consideradas en la Suma de Cuadrados del Error \[SSE=\sum{(Y_i-\hat{Y})^2}\]

Finalmente, la diferencia entre las sumas de cuadrados anteriores se toma en cuenta con las desviaciones del valor ajustado \(\hat{Y}\) con respecto al promedio de las observaciones \(\bar{Y}\) a través de la Suma de Cuadrados de la Regresión \[SSR=\sum{(\hat{Y}_i-\bar{Y})^2}\]

Notemos que \[(Y_i-\bar{Y})=(\hat{Y}_i-\bar{Y})+(Y_i-\hat{Y}_i)\]

Es interesante ver que esta descomposición también es válida con las sumas de cuadrados: \[SSTO=SSR+SSE\]

Veamos: \[ \begin{align*} SSTO&=\sum{(Y_i-\bar{Y})^2}\\ &=\sum{[(\hat{Y}_i-\bar{Y})+(Y_i-\hat{Y}_i)]^2}\\ &=\sum{[(\hat{Y}_i-\bar{Y})^2+(Y_i-\hat{Y}_i)^2+2(\hat{Y}_i-\bar{Y})(Y_i-\hat{Y}_i)]}\\ &=\sum{(\hat{Y}_i-\bar{Y})^2}+\sum{(Y_i-\hat{Y}_i)^2}\\ &=SSR+SSE \end{align*} \]

Notar que \[\sum{2(\hat{Y}_i-\bar{Y})(Y_i-\hat{Y}_i)}=2\sum{\hat{Y}_ie_i}-2\bar{Y}\sum{e_i}=0\]

Modelo ANOVA de un factor

El modelo ANOVA asume que:

  • La distribución de probabilidad para cada nivel del factor es Normal.
  • Para cada nivel del factor, la varianza es la misma.
  • Las respuestas para cada nivel son muestras aleatorias de su respectiva distribución, independientes de las respuestas en los otros niveles del factor.

El análisis de las distribuciones de probabilidad correspondientes a los niveles del factor se realiza, en general, en dos etapas:

  • Determinar si las medias de los niveles son iguales o no.
  • Si hay diferencias entre los niveles, examinar cómo difieren y qué implicancias tienen esas diferencias.

Notación: \[ \begin{align*} r&=\text{# de niveles del factor en estudio.}\\ n_i&=\text{# de casos para el nivel }i, i=1,2,\ldots,r\\ n_T&=\text{# de casos en el estudio }n_T=\sum{n_i}\\ Y_{ij}&=\text{Valor de la respuesta en el }j\text{-esimo ensayo del nivel }i \end{align*} \]

El modelo ANOVA se escribe así: \[ \begin{align*} Y_{ij}&=\mu_i+\varepsilon_{ij}\text{ con }\varepsilon_{ij}\text{ i.i.d. }N(0,\sigma^2)\\ \mu_i&=E(Y_{ij})\ j=1,\ldots,n_i;\ i=1,\ldots,r\\ Var(Y_{ij})&=\sigma^2 \end{align*} \]

Matricialmente, si tenemos \(r=3\) tratamientos con \(n_1=n_2=n_3=2\), entonces el modelo se puede representar por \(Y=X\beta+\varepsilon\) con

\[ Y=\left[\begin{array}{l}Y_{11}\\Y_{12}\\Y_{21}\\Y_{22}\\Y_{31}\\Y_{32}\end{array}\right] X=\left[\begin{array}{lll}1&0&0\\1&0&0\\0&1&0\\0&1&0\\0&0&1\\0&0&1\end{array}\right] \beta=\left[\begin{array}{l}\mu_1\\\mu_2\\\mu_3\end{array}\right] \varepsilon=\left[\begin{array}{l}\varepsilon_{11}\\\varepsilon_{12}\\\varepsilon_{21}\\\varepsilon_{22}\\\varepsilon_{31}\\\varepsilon_{32}\end{array}\right] \]

\[Var(Y)=Var(\varepsilon)=\sigma^2I_{n_T\times n_T}\]

  • En estudios observacionales, las medias \(\mu_i\) corresponden a las medias de las poblaciones correspondientes a cada nivel
  • En estudios experimentales, la media \(\mu_i\) representa la respuesta media que se obtiene al aplicar el tratamiento \(i\) a todas las unidades en la población de unidades experimentales acerca de la cual queremos hacer {}.
  • Es decir, la diferencia es que en el segundo caso podemos atribuir la diferencia al tratamiento respectivo, mientras que en el primer caso la diferencia tiene sólo un valor descriptivo.
  • Denotaremos por \(Y_{i.}\) a la suma total de observaciones para el nivel \(i\)

\[Y_{i.}=\sum_{j=1}^{n_i}{Y_{ij}}\] \(\bar{Y}_{i.}\) es la media muestral para el nivel \(i\) \[\bar{Y}_{i.}=\frac{Y_{i.}}{n_i}\]

\(Y_{..}\) es la suma total para todas las observaciones del estudio \[Y_{..}=\sum_{i=1}^{r}\sum_{j=1}^{n_i}{Y_{ij}}\]

La media general de todas las observaciones es \[\bar{Y}_{..}=\frac{Y_{..}}{n_T}=\sum_{i=1}^{r}{\frac{n_i}{n_T}\bar{Y}_i}\]

Para mínimos cuadrados, la cantidad a minimizar es

\[ \begin{align*} Q&=\sum_i\sum_j{(Y_{ij}-\mu_i)^2}\\ &=\sum_j{(Y_{1j}-\mu_1)^2}+\sum_j{(Y_{2j}-\mu_2)^2}+\ldots+\sum_j{(Y_{rj}-\mu_r)^2} \end{align*} \]

O sea que \(Q\) se puede minimizar a través de cada componente por separado, lo que a su vez se realiza con \(\hat{\mu}_i=\bar{Y}_{i.}\). Es decir,

\[\hat{Y}_{ij}=\bar{Y}_{i.}\]

Los mismos estimadores son obtenidos al maximizar la verosimilitud

\[L(\mu_1,\ldots,\mu_r,\sigma^2)=\frac{1}{(2\pi\sigma^2)^{n/2}}\exp\left[-\frac{1}{2\sigma^2}\sum_i\sum_j{(Y_{ij}-\mu_i)^2}\right]\]

Los residuos están definidos por

\[e_{ij}=Y_{ij}-\hat{Y}_{ij}=Y_{ij}-\bar{Y}_i\] y por lo tanto

\[\sum_j{e_{ij}=0}\ i=1,2,\ldots,r\]

Suma de cuadrados total

\[SSTO=\sum_i\sum_j{(Y_{ij}-\bar{Y}_{..})^2}\]

Suma de cuadrados de los tratamientos (variabilidad entre tratamientos)

\[SSTR=\sum_i{n_i(\bar{Y}_{i.}-\bar{Y}_{..})^2}\]

Suma de cuadrados del error (variabilidad dentro de cada tratamiento)

\[SSE=\sum_i\sum_j{(Y_{ij}-\bar{Y}_{i.})^2}\]

Y tenemos \[SSTO=SSTR+SSE\]

La SSTO tiene \(n_T-1\) grados de libertad, ya que hay \(n_T\) desviaciones \((Y_{ij}-\bar{Y}_{..})\) y se pierde un grado de libertad porque

\[\sum\sum{(Y_{ij}-\bar{Y}_{..})}=0\]

La SSTR tiene \(r-1\) grados de libertad, ya que hay \(r\) desviaciones \((\bar{Y}_{i.}-\bar{Y}_{..})\) con

\[\sum{n_i(\bar{Y}_{i.}-\bar{Y}_{..})}=0\]

La SSE tiene \(n_T-r\) grados de libertad asociados. Esto es porque para cada nivel \(i\) el término \(\sum_j{(Y_{ij}-\bar{Y}_{i.})^2}\) tiene \(n_i-1\) g.l. y luego al sumar \((n_1-1)+\ldots+(n_r-1)\) obtenemos \(n_T-r\)

Tabla ANOVA

Tabla ANOVA

Trabajo en Clase: Calcular “a mano” los coeficientes estimados del modelo ANOVA y la Tabla ANOVA

MIKTEX MACTEX

Test F para igualdad de medias

Queremos docimar

\[ \begin{align*} H_0:&\mu_1=\mu_2=\ldots=\mu_r=\mu_0\\ H_1:&\text{No todos los }\mu_i\text{ son iguales.} \end{align*} \]

Bajo \(H_0\), \[F^*=\frac{MSTR}{MSE}\sim F_{r-1,n_T-r}\]

Si \(r=2\), entonces el test \(F^*\) con \(1\) y \(n_T-2\) grados de libertad es equivalente al test T de dos colas \[t^*=\frac{\bar{Y}_1-\bar{Y}_2}{\sqrt{S^2\left(\frac{1}{n_1}+\frac{1}{n_2}\right)}}\] con \[S^2=\frac{(n_1-1)S^2_1+(n_2-1)S^2_2}{n_1+n_2-2}\]

Mediante el test \(t^*\) se puede hacer pruebas de una cola o dos colas.

Si denotamos por (F) al modelo completo y (R) una versión reducida, entonces el test lineal general es

\[F^*=\frac{SSE(R)-SSE(F)}{g.l.(R)-g.l.(F)}\big/\frac{SSE(F)}{g.l.(F)}\]

En este caso particular, el modelo reducido es \[Y_{ij}=\mu_0+\varepsilon_{ij}\] y \[SSE(R)=\sum\sum{(Y_{ij}-\hat{Y}_{ij})^2}=\sum\sum{(Y_{ij}-\bar{Y}_{..})^2}=SSTO\]

y luego \[F^*=\frac{SSTO-SSE}{(n_T-1)-(n_T-r)}\big/\frac{SSE}{n_T-r}=\frac{SSTR}{r-1}\big/\frac{SSE}{n_T-r}=\frac{MSTR}{MSE}\]

Parametrización alternativa

A veces interesa representar las diferencias entre tratamientos en términos de un nivel de comparación \(\mu_0\). Es decir \[Y_{ij}=\mu_0+\tau_i+\varepsilon_{ij}\]

con \[\tau_i=\mu_i-\mu_0\]

Por ejemplo \[\mu_0=\sum_{i=1}^r{\omega_i\mu_i}\] con \[\sum_{i=1}^r{\omega_i}=1\]

Los \(\tau_i\) quedan sujetos a la restricción \[\sum{\omega_i\tau_i}=0\]

Entonces \(H_0:\mu_1=\mu_2=\ldots=\mu_r=\mu_0\) es equivalente a \[H_0:\tau_1=\tau_2=\ldots=\tau_r=0\] y el test \(F^*\) es el mismo.

Supongamos \[\omega_i=\frac{n_i}{n_T}\]

Entonces \(\sum{\frac{n_i}{n_T}}=0\) y podemos definir \[\tau_r=-\frac{n_1}{n_r}\tau_1-\frac{n_2}{n_r}\tau_2-\ldots-\frac{n_{r-1}}{n_r}\tau_{r-1}\]

y luego \[Y_{ij}=\mu_0+\tau_1X_{ij1}+\tau_2X_{ij2}+\ldots+\tau_{r-1}X_{ij(r-1)}+\varepsilon_{ij}\] con \[ \begin{align*} X_{ij1}&=\left\{\begin{array}{ll}1&\text{ nivel }1\\-\frac{n_1}{n_r}&\text{ nivel }r\\0&\text{ e.o.c.}\end{array}\right.\\ \cdots X_{ij(r-1)}&=\left\{\begin{array}{ll}1&\text{ nivel }r-1\\-\frac{n_{r-1}}{n_r}&\text{ nivel }r\\0&\text{ e.o.c.}\end{array}\right. \end{align*} \]

Supongamos que tenemos \(r=3\) niveles con \(n_1=n_2=n_3=2\) repeticiones. Entonces, matricialmente, el modelo es \[Y=X\beta+\varepsilon\] con \[ Y=\left[\begin{array}{c}Y_{11}\\Y_{12}\\Y_{21}\\Y_{22}\\Y_{31}\\Y_{32}\end{array}\right] X=\left[\begin{array}{ccc}1&1&0\\1&1&0\\1&0&1\\1&0&1\\1&-\frac{n_1}{n_3}&-\frac{n_2}{n_3}\\1&-\frac{n_1}{n_3}&-\frac{n_1}{n_3}\end{array}\right] \beta=\left[\begin{array}{c}\mu_0\\\tau_1\\\tau_2 \end{array}\right] \varepsilon=\left[\begin{array}{c}\varepsilon_{11}\\\varepsilon_{12}\\\varepsilon_{21}\\\varepsilon_{22}\\\varepsilon_{31}\\\varepsilon_{32}\end{array}\right] \]

Otras parametrizaciones en R

Celda de referencia

Referencia en el tratamiento 1 (lo que ya teníamos)

summary(modelo1)
## 
## Call:
## lm(formula = Venta ~ Diseño, data = datos)
## 
## 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)   14.600      1.452  10.053 4.66e-08 ***
## Diseño2       -1.200      2.054  -0.584   0.5677    
## Diseño3        4.900      2.179   2.249   0.0399 *  
## Diseño4       12.600      2.054   6.135 1.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.248 on 15 degrees of freedom
## Multiple R-squared:  0.7881, Adjusted R-squared:  0.7457 
## F-statistic: 18.59 on 3 and 15 DF,  p-value: 2.585e-05
stats::anova(modelo1)
## Analysis of Variance Table
## 
## Response: Venta
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Diseño     3 588.22 196.074  18.591 2.585e-05 ***
## Residuals 15 158.20  10.547                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Podemos pedir los valores ajustados del modelo

predict(object = modelo1, newdata = tibble(Diseño = factor(1 : 4)), type = "response"
)
##    1    2    3    4 
## 14.6 13.4 19.5 27.2

Veamos qué pasa si ponemos como referencia el tratamiento 2

contrasts(datos$Diseño) <- contr.treatment(nlevels(datos$Diseño), base = 1)
contrasts(datos$Diseño) #Entrega la matriz de CONTRASTES (codificación de variables indicadoras)
##   2 3 4
## 1 0 0 0
## 2 1 0 0
## 3 0 1 0
## 4 0 0 1
contrasts(datos$Diseño) <- contr.treatment(nlevels(datos$Diseño), base = 2)
contrasts(datos$Diseño) #Esta vez la matriz de contrastes tiene como referencia el tratamiento 2
##   1 3 4
## 1 1 0 0
## 2 0 0 0
## 3 0 1 0
## 4 0 0 1
modelo1r2 <- lm(Venta ~ Diseño, data = datos)
summary(modelo1r2)
## 
## Call:
## lm(formula = Venta ~ Diseño, data = datos)
## 
## 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)   13.400      1.452   9.226 1.43e-07 ***
## Diseño1        1.200      2.054   0.584   0.5677    
## Diseño3        6.100      2.179   2.800   0.0135 *  
## Diseño4       13.800      2.054   6.719 6.88e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.248 on 15 degrees of freedom
## Multiple R-squared:  0.7881, Adjusted R-squared:  0.7457 
## F-statistic: 18.59 on 3 and 15 DF,  p-value: 2.585e-05
stats::anova(modelo1r2)
## Analysis of Variance Table
## 
## Response: Venta
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Diseño     3 588.22 196.074  18.591 2.585e-05 ***
## Residuals 15 158.20  10.547                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict(object = modelo1r2, newdata = tibble(Diseño = factor(1 : 4)), type = "response"
)
##    1    2    3    4 
## 14.6 13.4 19.5 27.2

Ahora veamos la parametrización para usar como referencia la media general

contrasts(datos$Diseño) <- contr.sum(nlevels(datos$Diseño))
contrasts(datos$Diseño) #Esta vez la matriz de contrastes está definida por la restricción de que la suma de los coeficientes sea 0
##   [,1] [,2] [,3]
## 1    1    0    0
## 2    0    1    0
## 3    0    0    1
## 4   -1   -1   -1
modelo1rs <- lm(Venta ~ Diseño, data = datos)
summary(modelo1rs)
## 
## Call:
## lm(formula = Venta ~ Diseño, data = datos)
## 
## 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)  18.6750     0.7485  24.949 1.25e-13 ***
## Diseño1      -4.0750     1.2708  -3.207 0.005884 ** 
## Diseño2      -5.2750     1.2708  -4.151 0.000854 ***
## Diseño3       0.8250     1.3706   0.602 0.556221    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.248 on 15 degrees of freedom
## Multiple R-squared:  0.7881, Adjusted R-squared:  0.7457 
## F-statistic: 18.59 on 3 and 15 DF,  p-value: 2.585e-05
mean(datos$Venta)
## [1] 18.63158
stats::anova(modelo1rs)
## Analysis of Variance Table
## 
## Response: Venta
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Diseño     3 588.22 196.074  18.591 2.585e-05 ***
## Residuals 15 158.20  10.547                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict(object = modelo1rs, newdata = tibble(Diseño = factor(1 : 4)), type = "response"
)
##    1    2    3    4 
## 14.6 13.4 19.5 27.2

Notar que la referencia no es exactamente la media general ya que el diseño es desbalanceado. En realidad es la media de las medias

mean(predict(object = modelo1rs, newdata = tibble(Diseño = factor(1 : 4)), type = "response"
)
)
## [1] 18.675

Prueba de homogeneidad de varianzas

# library(car)
datos %>% 
  leveneTest(Venta ~ Diseño, data = .)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.2417 0.8659
##       15

Contrastes

Un contraste se define como una combinación lineal de las medias \(\mu_i\) de los niveles del factor, en que los coeficientes suman cero.

\[L=\sum_{i=1}^r{c_i\mu_i}\ \text{con }\sum_{i=1}^r{c_i=0}\]

Un estimador insesgado de un contraste \(L\) es \[\hat{L}= \sum_{i=1}^r{c_i\bar{Y}_{i.}}\]

con varianza \(Var(\hat{L})=\sigma^2\sum_{i=1}^r{\frac{c_i^2}{n_i}}\). Un estimador insesgado de esta varianza es

\[S^2(\hat{L})=MSE\sum_{i=1}^r{\frac{c_i^2}{n_i}}\]

Luego

\[\frac{\hat{L}-L}{S(\hat{L})}\]

tiene distribución \(t\) con \(n_T-r\) grados de libertad. Basado en esto se realiza el test de hipótesis

\[ \begin{eqnarray*} H_0:L&=&0\\ H_1:L&\neq&0 \end{eqnarray*} \]

Para combinaciones lineales de las medias, se aplica lo mismo.

Los procedimientos para comparar medias entre tratamientos que hemos visto hasta ahora tienen dos limitaciones importantes:

  • El nivel de confianza \(1-\alpha\) se aplica sólo a una estimación particular, no a una serie de estimaciones. Lo mismo para el error de tipo I, \(\alpha\).

  • El nivel de confianza \(1-\alpha\) es apropiado sólo si el estimador o el test no es resultado de observar los datos.

La solución a estas limitaciones es usar un test de comparaciones múltiples, que incluya en forma conjunta todas las comparaciones de interés potencial.

Comparación de Tukey

En este procedimiento, el interés está puesto en todas las combinaciones posibles de a pares, es decir, las pruebas del tipo

\[ \begin{eqnarray*} H_0:\mu_i-\mu_{i'}&=&0\\ H_1:\mu_i-\mu_{i'}&\neq&0 \end{eqnarray*} \]

Cuando todos los tamaños de muestra son iguales, el nivel de significancia global del método de Tukey es exactamente \(\alpha\). Este se reduce cuando los tamaños muestrales no son iguales, lo que lo convierte en una prueba conservadora.

Supongamos que tenemos \(r\) observaciones \(Y_1,\ldots,Y_r\) i.i.d. \(N(\mu,\sigma^2)\). Sea

\[\omega=\max_i(Y_i)-\min_i{Y_i}.\]

Sea \(S^2\) un estimador de la varianza \(\sigma^2\) basado en \(\nu\) grados de libertad e independiente de los \(Y_i\) (en nuestro caso \(\sigma^2=MSE\) y \(\nu=n_T-r\)). Entonces el cuociente

\[q(r,\nu)=\frac{\omega}{S}\]

se denomina rango studentizado, y su distribución ha sido tabulada.

Se construye intervalos de confianza múltiples para todas las comparaciones de a pares \(D=\mu_i-\mu_{i'}\) con nivel de confianza global de al menos \(1-\alpha\) basado en

\[\hat{D}\pm T\ S(\hat{D})\]

con

\[ \begin{eqnarray*} \hat{D}&=&\bar{Y}_{i.}-\bar{Y}_{i'.}\\ S^2(\hat{D})&=&MSE\left(\frac{1}{n_i}+\frac{1}{n_{i'}}\right)\\ T&=&\frac{1}{\sqrt{2}}q(1-\alpha;r,n_T-r) \end{eqnarray*} \]

Equivalentemente, se puede efectuar los tests de hipótesis correspondientes mediante el estadístico de prueba

\[q^*=\frac{\sqrt{2}\hat{D}}{S(\hat{D})}\]

el cual se compara con el percentil \(q(1-\alpha;r,n_T-r)\).

Comparación de Scheffé

Esta prueba se usa cuando el interés se centra en todos los posibles contrastes entre los niveles del factor, para pruebas de hipótesis del tipo

\[ \begin{eqnarray*} H_0:L&=&0\\ H_1:L&\neq&0 \end{eqnarray*} \]

El nivel de confianza global obtenido es exactamente \(1-\alpha\), sin importar si los tamaños muestrales en los niveles es igual o no. Los intervalos de confianza de Scheffé para la familia de contrastes \(L\) tienen la forma

\[\hat{L}\pm S\cdot S(\hat{L})\]

con

\[ \begin{eqnarray*} \hat{L}&=&\sum{c_i\bar{Y}_{i.}}\\ S^2(\hat{L})&=&MSE\sum\frac{c_i^2}{n_i}\text{ y}\\ S^2&=&(r-1)F(1-\alpha;r-1,n_T-r) \end{eqnarray*} \]

Para docimar directamente la hipótesis nula versus la alternativa, se ocupa el estadístico

\[F^*=\frac{\hat{L}^2}{(r-1)S^2(\hat{L})}\]

el cual se compara con el percentil \(F(1-\alpha;r-1,n_T-r)\).

Comparación de Bonferroni

Este método se aplica cuando la familia de interés es un conjunto de comparaciones de a pares, contrastes o combinaciones lineales en particular, especificada anteriormente al análisis estadístico.

Sea \(g\) el número de comparaciones de a pares, contrastes o combinaciones lineales de interés (en general, todas son combinaciones lineales). El método de Bonferroni asegura que los intervalos de confianza de todas las \(g\) combinaciones lineales \(L\) son correctos con un nivel de confianza de al menos \(1-\alpha\). Los intervalos de confianza están dados por

\[\hat{L}\pm t(1-\alpha/(2g);n_T-r)\cdot S(\hat{L}).\]

En forma equivalente, para la dócima de hipótesis

\[ \begin{eqnarray*} H_0:L&=&0\\ H_1:L&\neq&0 \end{eqnarray*} \]

se ocupa la estadistica de prueba

\[t^*=\frac{\hat{L}}{S(\hat{L})},\]

la cual se compara con \(t(1-\alpha/(2g);n_T-r)\) en un test de dos colas.

Comparación entre las pruebas de Tukey, Scheffé y Bonferroni

  • Si sólo interesa comparaciones de a pares (todas), el procedimiento de Tukey entrega intervalos de confianza más precisos que los de Scheffé y Bonferroni.

  • Si el test \(F\) para igualdad entre las medias de los niveles indica que las medias son distintas, la comparación múltiple de Scheffé encontrará al menos un contraste (de entre todos los posibles) que difiere significativamente de cero.

  • El procedimiento de Bonferroni es preferible al de Scheffé si el número de contrastes de interés es cercano o menor al número de niveles del factor.

  • En caso que dos o más procedimientos de comparación múltiple sean aplicables, se puede probar todos y elegir el más preciso, ya que la elección no depende de los datos observados.

Prueba de Tukey en R

modelo1.aov <- aov(Venta ~ Diseño, data = datos)
summary(modelo1.aov)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Diseño       3  588.2  196.07   18.59 2.58e-05 ***
## Residuals   15  158.2   10.55                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
stats::anova(modelo1.aov)
## Analysis of Variance Table
## 
## Response: Venta
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Diseño     3 588.22 196.074  18.591 2.585e-05 ***
## Residuals 15 158.20  10.547                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(modelo1.aov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Venta ~ Diseño, data = datos)
## 
## $Diseño
##     diff       lwr       upr     p adj
## 2-1 -1.2 -7.119758  4.719758 0.9352978
## 3-1  4.9 -1.378852 11.178852 0.1548895
## 4-1 12.6  6.680242 18.519758 0.0001013
## 3-2  6.1 -0.178852 12.378852 0.0582866
## 4-2 13.8  7.880242 19.719758 0.0000368
## 4-3  7.7  1.421148 13.978852 0.0142180

Contrastes personalizados

contr.diseño <- cbind(
  "2 - 1" = c(-1, 1, 0, 0),
  "3 - 1" = c(-1, 0, 1, 0),
  "4 - 1" = c(-1, 0, 0, 1),
  "3 - (1 + 2)" = c(-0.5, -0.5, 1, 0),
  "4 - (1 + 2)" = c(-0.5, -0.5, 0, 1)
)

comp.diseño <- glht(
  modelo1.aov,
  linfct = mcp(Diseño = t(contr.diseño))
)

summary(comp.diseño, test = univariate())
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = Venta ~ Diseño, data = datos)
## 
## Linear Hypotheses:
##                  Estimate Std. Error t value Pr(>|t|)    
## 2 - 1 == 0         -1.200      2.054  -0.584   0.5677    
## 3 - 1 == 0          4.900      2.179   2.249   0.0399 *  
## 4 - 1 == 0         12.600      2.054   6.135 1.91e-05 ***
## 3 - (1 + 2) == 0    5.500      1.921   2.863   0.0119 *  
## 4 - (1 + 2) == 0   13.200      1.779   7.421 2.15e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
summary(comp.diseño, test = adjusted("bonferroni"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = Venta ~ Diseño, data = datos)
## 
## Linear Hypotheses:
##                  Estimate Std. Error t value Pr(>|t|)    
## 2 - 1 == 0         -1.200      2.054  -0.584   1.0000    
## 3 - 1 == 0          4.900      2.179   2.249   0.1997    
## 4 - 1 == 0         12.600      2.054   6.135 9.55e-05 ***
## 3 - (1 + 2) == 0    5.500      1.921   2.863   0.0593 .  
## 4 - (1 + 2) == 0   13.200      1.779   7.421 1.07e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
datos_sim <- function(medias, d.e., k, n) {
  m = rep(rep(medias, each = n), k)
  ntrat = length(medias)
  tibble(Y = rnorm(n = n * k * ntrat, mean = m, sd = d.e.),
         Tratamiento = as.factor(rep(rep(1:3, each = n), k)),
         k = rep(1 : k, each = n * ntrat))
}

test_f <- function(y, x, alpha = 0.05) {
  s <- summary(
    lm(y ~ x)
  )
  tf <- s$fstatistic
  1 * ((1 - pf(tf[1], tf[2], tf[3])) < alpha)
}

potencia = function(medias, d.e., k, n, alpha = 0.05) {
  datos_sim(medias = medias, 
    d.e. = d.e., 
    k = k, 
    n = n) %>% 
    group_by(k) %>% 
    summarise(
      sig = test_f(Y, Tratamiento, alpha)
    ) %>% 
    summarise(Power = mean(sig)) %>% 
    mutate(              
      lci.95 = round(Power - 1.96 * sqrt(Power * (1 - Power) / k), 2),
      uci.95 = round(Power + 1.96 * sqrt(Power * (1 - Power) / k), 2),
      n = n,
      alpha = alpha)
}

Estudio métodos de extracción en hojas de ficus carica

Efecto de los métodos de extracción en las siguientes variables respuesta:

  • Polifenoles totales
  • Capacidad antioxidante
  • Variable por definir

Los tratamientos son:

  1. Convencional 1.1 con solvente Etanol 1.2 con solvente Metanol 1.3 con Agua (Control)
  2. Extracción Supercrítica 2.1 CO2 2.2 CO2 - Etanol 5% 2.3 CO2 - Etanol 7-10 %
  3. Ultrasonido 3.x Distintos niveles de Potencia (4), Tiempo (4) y Temperatura (3)

Estudio de potencia

Vamos a suponer ciertas diferencias en la variable respuesta de acuerdo a los distintos tratamientos, y vamos a evaluar empíricamente la potencia de una prueba ANOVA para diferentes tamaños muestrales (réplicas por tratamiento).

Sean \(Y_{ij}\) las observaciones de la variable respuesta para los tratamientos \(i = 1,2,\ldots, k\), con \(j = 1, 2, \ldots, n\) réplicas por tratamiento. Es decir, estamos suponiendo un diseño balanceado. Suponemos que

\[ Y_{ij} \sim N(\mu_i, \sigma^2) \]

En nuestro caso, vamos a considerar que \(Y\) representa los polifenoles totales luego de la extracción, y supongamos que tenemos \(k = 3\) tratamientos: 1.1, 1.2 y 1.3.

Simulación de datos

Vamos a suponer que la verdadera distribución de la variable respuesta está dada por \(\mu_1 = 4.545\), \(\mu_2 = 4.727\), \(\mu_3 = 6.909\), \(\sigma = 1\)

y vamos a ver qué capacidad tiene un modelo ANOVA de detectar las diferencias para distintos tamaños muestrales.

potencia(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 1000, n = 3, alpha = 0.05)
## # A tibble: 1 x 5
##   Power lci.95 uci.95     n alpha
##   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1   0.6  0.570   0.63     3  0.05
datos_sim(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 1000, n = 10) %>% 
  group_by(k) %>% 
  summarise(
      sig = test_f(Y, Tratamiento, 0.05)
    )
## # A tibble: 1,000 x 2
##        k   sig
##    <int> <dbl>
##  1     1     1
##  2     2     1
##  3     3     1
##  4     4     1
##  5     5     1
##  6     6     1
##  7     7     1
##  8     8     1
##  9     9     1
## 10    10     1
## # ... with 990 more rows
potencia(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 1000, n = 3, alpha = 0.01)
## # A tibble: 1 x 5
##   Power lci.95 uci.95     n alpha
##   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 0.256   0.23   0.28     3  0.01
potencia(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 1000, n = 5, alpha = 0.05)
## # A tibble: 1 x 5
##   Power lci.95 uci.95     n alpha
##   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 0.911   0.89   0.93     5  0.05
potencia(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 1000, n = 5, alpha = 0.01)
## # A tibble: 1 x 5
##   Power lci.95 uci.95     n alpha
##   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 0.706   0.68   0.73     5  0.01
potencia(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 1000, n = 10, alpha = 0.05)
## # A tibble: 1 x 5
##   Power lci.95 uci.95     n alpha
##   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1     1      1      1    10  0.05
potencia(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 1000, n = 10, alpha = 0.01)
## # A tibble: 1 x 5
##   Power lci.95 uci.95     n alpha
##   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 0.993   0.99      1    10  0.01
potencia(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 10000, n = 3, alpha = 0.05)
## # A tibble: 1 x 5
##   Power lci.95 uci.95     n alpha
##   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 0.600   0.59   0.61     3  0.05
potencia(medias = c(4.545, 4.727, 6.909), d.e. = 1, k = 1000, n = 3, alpha = 0.05)
## # A tibble: 1 x 5
##   Power lci.95 uci.95     n alpha
##   <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 0.614  0.580   0.64     3  0.05