Tema 8: Análisis de la varianza

class: center, middle, inverse, title-slide

.title[
# Tema 8: Análisis de la varianza
]
.date[
### <br><br><br>Estadística. Grado en Ciencias del Mar
]

---

## Planteamiento del problema: .blue[Ejemplo]

.pull-left[
Las condiciones de conservación del pescado se evalúan a través de la concentración de TVBN (_Total Volatile Base Nitrogen_). A mayor concentración de este elemento, peor es el estado de conservación de la pieza. Con objeto de determinar la temperatura que produce la menor concentración de TVBN, se eligen al azar 30 atunes recién pescados, todos de idéntico peso y caracterı́sticas generales. Se separan en tres grupos de 10 piezas cada uno. El primer grupo se congela a -4ºC, el segundo a -20ºC y el tercero a -40ºC. La tabla adjunta muestra la concentración de TVBN
en cada pieza después de 2 semanas de congelación. 
]

.pull-right[
<p style="padding: 2px 10px 2px 10px;"></p>
<table class="table" style="font-size: 18px; width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:right;font-weight: bold;color: black !important;background-color: darkgrey !important;"> Temperatura </th>
   <th style="text-align:right;font-weight: bold;color: black !important;background-color: darkgrey !important;"> -4ºC </th>
   <th style="text-align:right;font-weight: bold;color: black !important;background-color: darkgrey !important;"> -20ºC </th>
   <th style="text-align:right;font-weight: bold;color: black !important;background-color: darkgrey !important;"> -30ºC </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 18.3 </td>
   <td style="text-align:right;"> 11.7 </td>
   <td style="text-align:right;"> 16.64 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 15.92 </td>
   <td style="text-align:right;"> 12.87 </td>
   <td style="text-align:right;"> 17.83 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 18.71 </td>
   <td style="text-align:right;"> 11.77 </td>
   <td style="text-align:right;"> 19.01 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 17.92 </td>
   <td style="text-align:right;"> 12.23 </td>
   <td style="text-align:right;"> 17.33 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 15.66 </td>
   <td style="text-align:right;"> 13.62 </td>
   <td style="text-align:right;"> 17.06 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 17.14 </td>
   <td style="text-align:right;"> 13.24 </td>
   <td style="text-align:right;"> 18.04 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 15.21 </td>
   <td style="text-align:right;"> 14.02 </td>
   <td style="text-align:right;"> 17.51 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 19.92 </td>
   <td style="text-align:right;"> 13.66 </td>
   <td style="text-align:right;"> 19.11 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 17.61 </td>
   <td style="text-align:right;"> 12.27 </td>
   <td style="text-align:right;"> 17.75 </td>
  </tr>
  <tr>
   <td style="text-align:right;">  </td>
   <td style="text-align:right;"> 13.43 </td>
   <td style="text-align:right;"> 12.45 </td>
   <td style="text-align:right;"> 19.36 </td>
  </tr>
  <tr>
   <td style="text-align:right;font-weight: bold;color: white !important;background-color: rgba(40, 80, 228, 255) !important;"> Media (Sd) </td>
   <td style="text-align:right;font-weight: bold;color: white !important;background-color: rgba(40, 80, 228, 255) !important;"> 17.96 (0.92) </td>
   <td style="text-align:right;font-weight: bold;color: white !important;background-color: rgba(40, 80, 228, 255) !important;"> 12.78 (0.82) </td>
   <td style="text-align:right;font-weight: bold;color: white !important;background-color: rgba(40, 80, 228, 255) !important;"> 16.98 (1.92) </td>
  </tr>
</tbody>
</table>

]

---

## Planteamiento del problema: .blue[Ejemplo]

.pull-left[

Podemos representar estos datos mediante un boxplot:

![](tema8_Análisis_de_la_varianza_files/figure-html/unnamed-chunk-2-1.png)

]

.pull-right[

Obviamente las medias muestrales son distintas:

<table class="table" style="font-size: 18px; width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;">  </th>
   <th style="text-align:right;"> -4ºC </th>
   <th style="text-align:right;"> -20ºC </th>
   <th style="text-align:right;"> -30ºC </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Media (sd) </td>
   <td style="text-align:right;"> 17.96 (0.92) </td>
   <td style="text-align:right;"> 12.78 (0.82) </td>
   <td style="text-align:right;"> 16.98 (1.92) </td>
  </tr>
</tbody>
</table>

<br>

.resalta[
Pero: __¿Podemos deducir de estos datos que existe evidencia suficiente de que las medias poblacionales difieren entre sí?__
]

]

---

## Planteamiento del problema: .blue[Ejemplo]

En definitiva, nos planteamos el contraste de hipótesis:

<br>

.resalta[

`$$\large\begin{cases} 
H_{0}: & \mu_{1}=\mu_{2}=\mu_{3}\\
H_{1:} & \exists i,j:\,\mu_{i}\ne\mu_{j}
\end{cases}$$`

]

<br>

siendo `$\mu_1$`, `$\mu_2$` y `$\mu_3$` las respectivas medias poblacionales de TVBN a cada una de las tres temperaturas.

---

## Análisis de la varianza.

En el caso más elemental, se denomina .blue[ __análisis de la varianza__] _con un .blue[factor de variación]_ al contraste:

.resalta[

`$$\large\begin{cases} 
H_{0}: & \mu_{1}=\mu_{2}=\dots = \mu_{p}\\
H_{1:} & \exists i,j:\,\mu_{i}\ne\mu_{j}
\end{cases}$$`

]

<br>

* El .blue[ __factor de variación__] es la variable que define las poblaciones que se comparan. El factor de variación es siempre una variable categórica, y sus valores se denominan .blue[ __niveles__].

* En nuestro ejemplo el factor de variación es la temperatura, que toma tres niveles (-4ºC, -20ºC y -40ºC)

* Nótese que aunque el contraste se denomina _análisis de la varianza_ es realmente un .red[contraste de comparación de medias]. Pronto veremos la razón de este nombre.

---

## Análisis de la varianza.

* ¿Por qué no realizar varias comparaciones mediante el test de la t de Student que ya hemos visto?

* Imaginemos que `$p$` (el número de medias a comparar) es `$p=7$`. Si queremos comparar todas las medias dos a dos tendríamos que hacer `$\binom{7}{2}=21$` contrastes.

* Si cada contraste se realiza con un nivel de significación `$\alpha=5\%$`, ello significa que de cada 100 contrastes en los que `$H_0$` fuese cierta, dicha hipótesis se rechazaría en 5; o lo que es lo mismo, se rechazaría incorrectamente en uno de cada 20 contrastes.

* Por tanto, si `$H_0: \mu_1=\mu_2=\dots=\mu_p$` es cierta, al hacer los 21 tests para comparar estas 7 medias dos a dos, .blue[ __¡podemos estar casi seguros de que en alguno de los contrastes decidiríamos que las dos medias implicadas son distintas!__], lo que nos obligaría a rechazar `$H_0$`.

* Dicho de otra forma, la probabilidad de rechazar `$H_0$` siendo cierta se "infla" cuando se realizan muchos contrastes de la t de Student para comparar las distintas medias dos a dos.

---

## Análisis de la varianza: Procedimiento

Los datos disponibles para el análisis de la varianza normalmente son de la forma:

`$$\begin{array}{cccccc}
\underline{\textrm{Grupo 1}} & \underline{\textrm{Grupo 2}} &  & \underline{\textrm{Grupo }i} &  & \underline{\textrm{Grupo }p}\\
y_{11} & y_{21} & \ldots & y_{i1} & \ldots & y_{p1}\\
y_{11} & y_{21} & \ldots & y_{i1} & \ldots & y_{p1}\\
\vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\
y_{1n_{1}} & y_{2n_{2}} & \ldots & y_{in_{i}} & \ldots & y_{pn_{p}}
\end{array}$$`

<br>

* Tenemos `$p$` grupos

* En el grupo `$i$` hay `$n_i$` observaciones

* Llamaremos `$\bar{y}_i$` a la media del grupo `$i$`, e `$\bar{y}$` a la media de todos los valores observados.

---

## Análisis de la varianza: Procedimiento

Sean:

* `${N=\sum\limits _{i=1}^{p}n_{i}}$` el número total de observaciones.

* `${\large S_{i}^{2}=\frac{\sum\limits _{j=1}^{n_{i}}{(y_{ij}-\bar{y}_{i})^{2}}}{n_{i}-1}}$` la varianza de los valores dentro del grupo `$i$`

* `${\large S_{R}^{2}=\frac{{\sum\limits _{i=1}^{p}{\sum\limits _{j=1}^{n_{i}}{(y_{ij}-\bar{y}_{i})^{2}}}}}{{N-p}}=\frac{{\sum\limits _{i=1}^{p}{\left({n_{i}-1}\right)S_{i}^{2}}}}{{N-p}}}$` la .blue[ __Variabilidad dentro de los grupos__]; esta medida es un promedio de las varianzas dentro de cada grupo.

* `${\large S_{E}^{2}=\frac{{\sum\limits _{i=1}^{p}{n_{i}(\bar{y}_{i}-\bar{y})^{2}}}}{{p-1}}}$` la .blue[ __variabilidad entre grupos__]; esta cantidad mide la variabilidad entre las medias de los grupos.

---

## Análisis de la varianza: Procedimiento

Sea `$\sigma^2_i$` la varianza poblacional en el grupo `$i$`. Supondremos que se dan las siguientes condiciones:

* .blue[ __Homoscedasticidad:__] Todos los grupos tienen la misma variabilidad: `${ \sigma^2_1=\sigma^2_2=\dots=\sigma^2_p=\sigma^2}$`

* .blue[ __Normalidad:__] Dentro de cada grupo, la variable respuesta observada sigue una distribución normal: `$Y_i\approx N\left(\mu_i,\sigma\right)$`

Si se dan estas condiciones se puede probar que:

`$$F_{exp}=\frac{{S_{E}^{2}}}{{S_{R}^{2}}}=\frac{{\frac{1}{{p-1}}\sum\limits _{i=1}^{p}{n_{i}(\bar{Y}_{i}-\bar{Y})^{2}}}}{{\frac{1}{{N-p}}\sum\limits _{i=1}^{p}{\sum\limits _{j=1}^{n_{i}}{(Y_{ij}-\bar{Y}_{i})^{2}}}}}\approx F_{p-1,N-p}$$`
--

<p style="background-color: #ffff80;">Si `$H_0$` es cierta la variabilidad entre grupos `$S^2_E$` debe ser menor que la variabilidad dentro de los grupos `$S^2_R$`, y por tanto `$F_{exp}$` debe ser un valor pequeño.</p>

---

## Análisis de la varianza: Procedimiento

Por tanto la regla de decisión consistirá en determinar cuál es el mayor valor que puede tomar `$F_{exp}$` por puro azar y rechazar `$H_0$` si `$F_{exp}$` es mayor que dicho valor. Concretamente, si fijamos `$\alpha$` como nivel de significación, la regla de decisión es:

<br>

.resalta[

* Si `$F_{exp}\le F_{p-1,N-p,\alpha}\Longrightarrow$` Aceptar `$H_{0}:\mu_{1}=\mu_{2}=...=\mu_{p}$`

* Si `$F_{exp}>F_{p-1,N-p,\alpha}\Longrightarrow$` Rechazar `$H_{0}$` y concluir que `$\exists\;\;i,j:\mu_{i}\ne\mu_{j}$`

]

---

## Análisis de la varianza: .blue[Ejemplo]

En nuestro ejemplo de la conservación de atunes se tiene que:

<table class="table" style="font-size: 14px; width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;">  </th>
   <th style="text-align:right;">  </th>
   <th style="text-align:right;">  </th>
   <th style="text-align:right;">  </th>
   <th style="text-align:right;">  </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Temp </td>
   <td style="text-align:right;"> -4ºC </td>
   <td style="text-align:right;"> -20ºC </td>
   <td style="text-align:right;"> -40ºC </td>
   <td style="text-align:right;"> GLOBAL </td>
  </tr>
  <tr>
   <td style="text-align:left;"> media </td>
   <td style="text-align:right;"> 17.964 </td>
   <td style="text-align:right;"> 12.783 </td>
   <td style="text-align:right;"> 16.982 </td>
   <td style="text-align:right;"> 15.91 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> sd </td>
   <td style="text-align:right;"> 0.919 </td>
   <td style="text-align:right;"> 0.823 </td>
   <td style="text-align:right;"> 1.923 </td>
   <td style="text-align:right;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> var </td>
   <td style="text-align:right;"> 0.845 </td>
   <td style="text-align:right;"> 0.677 </td>
   <td style="text-align:right;"> 3.698 </td>
   <td style="text-align:right;">  </td>
  </tr>
</tbody>
</table>

Entonces:

* `$S_{E}^{2}=\frac{{\sum\limits _{i=1}^{p}{n_{i}(\bar{Y}_{i}-\bar{Y})^{2}}}}{{p-1}}=\frac{10\left(17.964-15.91\right)^{2}+10\left(12.783-15.91\right)^{2}+10\left(16.982-15.91\right)^{2}}{3-1}=75.73$`

* `$S_{R}^{2}=\frac{{\sum\limits _{i=1}^{p}{\left({n_{i}-1}\right)S_{i}^{2}}}}{{N-p}}=\frac{9\cdot0.845+9\cdot0.677+9\cdot3.698}{30-3}=1.74$`

`$$\left.\begin{array}{cc}
F_{exp}= & \frac{S_{E}^{2}}{S_{R}^{2}}=\frac{75.73}{1.74}=43.52\\
F_{p-1,N-p,\alpha}= & F_{2,27,0.05}\,=\,3.354
\end{array}\right\} \Rightarrow F_{exp}>F_{2,27,0.05}\Rightarrow\textrm{Se rechaza }H_{0}$$`

---

## Análisis de la varianza con R: .blue[Ejemplo]

Para realizar el análisis de la varianza con R, los datos deben guardarse en dos columnas, una con la temperatura y otra con el valor de TVBN:

.left-column[
![](imagenes/archivoNitrogeno.png)
]

.right-column[

El código en R para el análisis de estos datos es:

```r
library(readxl)
atunes <- read_excel("datos/nitrogenoVolatil.xlsx")
summary(aov(TVBN~factor(Temp), data=atunes))
```

```
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(Temp)  2 151.46   75.73   43.52 3.58e-09 ***
## Residuals    27  46.99    1.74                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

]

Como el p-valor (3.58e-09) es menor que el nivel de significación 0.05 concluimos que existe evidencia suficiente para rechazar `$H_0$`

</p>

---

## Análisis de la varianza: Presentación de resultados

Normalmente, los resultados de un análisis de la varianza se presentan en una tabla análoga a la que ha mostrado R:

<br>

`$${\small\begin{array}{cccccc}
\underline{\textrm{Fuente de variación}} & \underline{\textrm{g.l.}} & \underline{\textrm{Suma de cuadrados}} & \underline{\textrm{Cuadrados medios}} & \underline{F} & \underline{\textrm{P-valor}}\\
\textrm{Factor} & p-1 & \sum\limits _{i=1}^{p}n_{i}(\bar{y}_{i}-\bar{y})^{2} & S_{E}^{2}=\frac{1}{p-1}\sum\limits _{i=1}^{p}n_{i}(\bar{Y}_{i}-\bar{Y})^{2} & S_{E}^{2}/S_{R}^{2} & p\\
\textrm{Residual} & N-p & \sum\limits _{i=1}^{p}\sum\limits _{j=1}^{n_{i}}(y-\bar{y}_{i})^{2} & S_{R}^{2}=\frac{1}{N-p}\sum\limits _{i=1}^{p}\sum\limits _{j=1}^{n_{i}}(Y_{ij}-\bar{Y}_{i})^{2}
\end{array}}$$`

<br>

El p-valor se calcula como:

`$$p=P\left(F_{p-1,N-p}>\frac{S^2_E}{S^2_R}\right)$$`

---

## Análisis de la varianza: Validación

Para validar el resultado de un análisis de la varianza hay que comprobar que se cumplen las condiciones en que se basa:

* .blue[Homoscedasticidad:] se comprueba mediante el test de Levene, cuya hipótesis nula es que los datos son homoscedásticos.

* .blue[Normalidad:] se comprueba mediante el test de Shapiro-Wilk, cuya hipótesis nula es que los datos son normales.

Ambos contrastes son laboriosos de realizar. Para llevarlos a cabo utilizaremos R y tomaremos la decisión basándonos en el p-valor de cada contraste:

* Si `$p-valor<\alpha$` se rechaza `$H_0$`

* Si `$p-valor\ge \alpha$`, se acepta `$H_0$`

Veamos a continuación la aplicación de estos contrastes en nuestro ejemplo.

---

## Análisis de la varianza: Validación con R

* __Homoscedasticidad:__

```r
library(car)
leveneTest(TVBN~factor(Temp), data=atunes)
```

```
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)  
## group  2  3.2265 0.0554 .
##       27                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

--
<br>

.resalta[
Como el p-valor es mayor que 0.05 podemos aceptar la hipótesis de homoscedasticidad.
]

---

## Análisis de la varianza: Validación con R

* __ Normalidad:__

.pull-left[

+ Contraste gráfico

```r
modelo=aov(TVBN~factor(Temp),data=atunes)
qqPlot(residuals(modelo))
```

![](tema8_Análisis_de_la_varianza_files/figure-html/unnamed-chunk-8-1.png)

```
## [1] 10  8
```
]

.pull-right[

+ Contraste numérico

```r
shapiro.test(residuals(modelo))
```

```
## 
## 	Shapiro-Wilk normality test
## 
## data:  residuals(modelo)
## W = 0.9783, p-value = 0.7788
```
  
  .resalta[
  Como el p-valor es mayor que 0.05 podemos aceptar la hipótesis de normalidad
  ]
]

---

## Análisis de la varianza: Contrastes a posteriori

Si en el contraste de la F del análisis de la varianza se acepta la existencia de diferencias significativas entre algunas de las medias, resulta de interés determinar qué poblaciones tienen medias diferentes y cuál es la magnitud de dichas diferencias. Para contrastar si las medias de las poblaciones r y s son iguales ó distintas:

`$$\begin{cases}
H_{0}: & \mu_{r}-\mu_{s}=0\\
H_{1}: & \mu_{r}-\mu_{s}\ne0
\end{cases}$$`

Puede utilizarse el .blue[ __Test de Scheffe__]

`$$\textrm{Si }\left|\frac{\bar{y}_{r}-\bar{y}_{s}}{S_{R}\sqrt{(p-1)\cdot\left({\frac{1}{{n_{r}}}+\frac{1}{{n_{s}}}}\right)}}\right|\le\sqrt{F_{p-1,N-p,\alpha}}	\Longrightarrow\textrm{Aceptar  } H_{0}$$`
`$$\textrm{En caso contrario}	\Longrightarrow\textrm{Rechazar } H_{0}$$`

---

## Análisis de la varianza: Contrastes a posteriori

Además un intervalo de confianza para `$\mu_{r}-\mu_{s}$` es:

`$$\mu_{r}-\mu_{s}\in\left[{\bar{y}_{r}-\bar{y}_{s}\pm S_{R}\sqrt{(p-1)\cdot F_{p-1,N-p,\alpha}\cdot\left({\frac{1}{{n_{r}}}+\frac{1}{{n_{s}}}}\right)}}\right]$$`
---

## Análisis de la varianza: Test de Scheffe con R

```r
library(DescTools)
ScheffeTest(modelo)
```

```
## 
##   Posthoc multiple comparisons of means: Scheffe Test 
##     95% family-wise confidence level
## 
## $`factor(Temp)`
##           diff    lwr.ci     upr.ci    pval    
## -20--40 -4.881 -6.434821 -3.3271794 5.5e-08 ***
## -4--40  -0.982 -2.535821  0.5718206  0.2788    
## -4--20   3.899  2.345179  5.4528206 3.0e-06 ***
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## Violación de los supuestos del Análisis de la varianza

1. Si se cumple el supuesto de normalidad, pero los datos __no son homoscedásticos__, se puede usar el __test de Welch__. En este caso se define:

`$$F_{exp}=\frac{\frac{1}{p-1}\sum_{j=1}^{p}w_{j}\left(\overline{Y}_{j}-\overline{Y}_{w}\right)^{2}}{1+\frac{2\left(p-2\right)}{p^{2}-1}\sum_{j=1}^{p}\left(\frac{1}{n_{j}-1}\right)\left(1-\frac{w_{j}}{w}\right)^{2}}$$`
donde:

`$$w_{j}=\frac{n_{j}}{s_{j}^{2}}\,\,\,\,\,\,\,\,\,\,w=\sum_{j=1}^{p}w_{j}\,\,\,\,\,\,\,\,\,\,\overline{Y}_{w}=\frac{\sum_{j=1}^{p}w_{j}\overline{Y}_{j}}{w}$$`
Cuando `$H_0$` es cierta `$F_{exp}\approx F_{p-1,d}$` con `$d=\frac{p^{2}-1}{3\sum_{j=1}^{p}\left(\frac{1}{n_{j}-1}\right)\left(1-\frac{w_{j}}{w}\right)^{2}}$`. La regla de decisión es rechazar `$H_0$` si `$F_{exp}>F_{p-1,d,\alpha}$`

---

2. Si no se cumplen el supuesto de normalidad ni el de homoscedasticidad, se puede usar el __test de Kruskal-Wallis__. En este caso se define:

`$$H=\frac{12}{N\left(N+1\right)}\sum_{j=1}^{p}\frac{R_{j}^{2}}{n_{j}}-3\left(N+1\right)$$`
siendo `$N$` el número total de observaciones y `$R_j$` la suma de rangos dentro de cada grupo. Para calcular los `$R_j$` se ordenan todos los valores de `$Y$` (sin importar el grupo) de menor a mayor; el rango de cada observación es su número de orden.

Cuando `$H_0$` es cierta `$H\approx \chi^2_{p-1}$`. La regla de decisión es entonces rechazar `$H_0$` cuando `$H>\chi^2_{p-1,\alpha}$`

---

## Contrastes a posteriori cuando no se cumplen los supuestos de Análisis de la varianza

1. Si se utiliza el test de Welch se pueden utilizar las pruebas Post-Hoc de Games-Howell

<br>

2. Si se utiliza el test de Kruskal-Wallis se pueden usar las comparaciones múltiples de Dwass-Steel-CRitchlow-Fligner