C08.2 ANOVA (un solo factor).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![imagenes](logo.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anova de un factor\n",
    "\n",
    "Como hemos dicho, el ANOVA de un solo factor es una prueba de hipótesis para saber si las medias de cierta característica  $X$  de la población son iguales en cada una de las $k$ clases (niveles) en que se clasifica la población de acuerdo a una característica $F$ (factor).\n",
    "\n",
    "Los supuestos para poderlo aplicar son los siguientes:\n",
    "\n",
    "* Las muestras de cada subpoblación son independientes interna y externamente.\n",
    "* Cada una de las subpoblaciones son gaussianas.\n",
    "* Todas las subpoblaciones tienen la misma varianza común $\\sigma^2$ (**homocedasticidad**).\n",
    "\n",
    "De esta manera, **imagina** que tu tabla se ve de la siguiente forma:\n",
    "\n",
    "$F_1$|$F_2$|...|$F_k$\n",
    ":--:|:--:|:--:|:--:\n",
    "$X_1^{(1)}$|$X_1^{(2)}$|...|$X_1^{(k)}$\n",
    "$X_2^{(1)}$|$X_2^{(2)}$|...|$X_2^{(k)}$\n",
    "$X_3^{(1)}$|$X_3^{(2)}$|...|$X_3^{(k)}$\n",
    "$\\vdots$|$\\vdots$|$\\vdots$|$\\vdots$\n",
    "$X_{n_1}^{(1)}$|$X_{n_2}^{(2)}$|...|$X_{n_k}^{(k)}$\n",
    "\n",
    "Observa que el total de elementos de la muestra pertencientes a la primer subpoblación es $n_1$; el total de elementos de la muestra pertencientes a la segunda subpoblación es $n_2$;...;el total de elementos de la muestra pertencientes a la $k$-ésima subpoblación es $n_k$. **No necesariamente se tiene $n_1=n_2=...=n_k$** en general.\n",
    "\n",
    "De esta manera, $X_{j}^{(i)}$ es el valor de la característica $X$ correspondiente al individuo $j$ del nivel $i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se definen las siguientes variables aleatorias:\n",
    "\n",
    "Variabla aleatoria|Descripción|Significado\n",
    ":--:|:--|:--\n",
    "$T^{(i)}$|suma del grupo $i$| suma de la columna $i$\n",
    "$\\overline{X^{(i)}}$|media del grupo $i$| media de la columna $i$\n",
    "$T$|suma de las $T^{(i)}$|suma de todos los valores\n",
    "$\\overline{X}$|$T$/total de datos|media de todos los valores\n",
    "\n",
    "Si $\\mu$ es la media global de toda la población y $\\mu_i$ es la media poblacional del grupo $i$, entonces $$X_j^{(i)}-\\mu=(X_j^{(i)}-\\mu_i)+(\\mu_i-\\mu)$$\n",
    "\n",
    "Si tomamos $\\overline{X^{(i)}}$ como estimador de $\\mu_i$ y $\\overline{X}$ como estimador de $\\mu$, elevando al cuadrado la identidad anterior tenemos: \n",
    "\n",
    "**Identidad de la suma de cuadrados para ANOVA**\n",
    "\n",
    "$$\\underbrace{\\sum_{i=1}^k\\sum_{j=1}^{n_i}(X_j^{(i)}-\\overline{X})^2}_{SS_{Total}}=\\underbrace{\\sum_{i=1}^kn_i(\\overline{X^{(i)}}-\\overline{X})^2}_{SS_{Tra}}+\\underbrace{\\sum_{i=1}^k\\sum_{j=1}^{n_i}(X_j^{(i)}-\\overline{X^{(i)}})^2}_{SS_{Err}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estadístico de contraste\n",
    "\n",
    "Partiendo de la igualdad anterior, se puede demostrar que $$(n-1)\\sigma^2_p=E\\left[SS_{Tra}\\right]+(n-k)\\sigma^2,$$ de modo que el valor promedio de $SS_{Tra}$ no depende de las submuestras. \n",
    "    \n",
    "Por otra parte, tomemos $MS_{Tra}=\\frac{SS_{Tra}}{k-1}$ y $MS_{Err}=\\frac{SS_{Err}}{n-k}$. Entonces es posible demostrar (pero es complicado) que\n",
    "\n",
    "$$\\begin{array}{rcl}E[MS_{Tra}]&=&\\sigma^2+\\frac{1}{k-1}\\sum_{i=1}^kn_i(\\mu_i-\\mu)^2\\\\E[MS_{Err}]&=&\\sigma^2\\end{array}$$\n",
    "\n",
    "De aquí concluimos dos cosas importantes:\n",
    "\n",
    "* $MS_{Err}$ es un estimador de la varianza común\n",
    "* $\\frac{MS_{Tra}}{MS_{Err}}\\ge 1$\n",
    "\n",
    "Si $H_0$ es cierta, entonces $\\frac{MS_{Tra}}{MS_{Err}}=1$ y si es falsa, entonces $\\frac{MS_{Tra}}{MS_{Err}}>1$.\n",
    "\n",
    "Toma $F$ como $$F=\\frac{MS_{Tra}}{MS_{Err}}.$$\n",
    "\n",
    "Si $H_0$ es cierta, la distribución de $F$ es lo que se conoce como una **distribución de Fisher** con $k-1$ y $n-k$ grados de libertad. Este es el estadístico de contraste de un ANOVA de un factor.\n",
    "\n",
    "Por lo tanto, **rechazamos $H_0$ si $F$ es bastante mayor que 1**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pasos del ANOVA\n",
    "\n",
    "Ahora, recordemos que para aplicar un ANOVA de un factor debe cumplirse que cada subpoblación es gaussiana y todas tienen la misma varianza. Además, una vez que tenemos el resultado, debemos conocer cuáles subplobaciones tienen varianzas diferentes. \n",
    "\n",
    "En **R**, para poder realizar ANOVA de forma automática se requiere, en primer lugar, que nuestra tabla se presente en formato largo: Es decir, de la forma\n",
    "\n",
    "Característica $X$|Factor $F$\n",
    ":--:|:--:\n",
    "$X_{1}^{(1)}$|1\n",
    "$X_{2}^{(1)}$|1\n",
    "$\\vdots$|$\\vdots$\n",
    "$X_{n_1}^{(1)}$|1\n",
    "$X_{1}^{(2)}$|2\n",
    "$X_{2}^{(2)}$|2\n",
    "$\\vdots$|$\\vdots$\n",
    "$X_{n_2}^{(2)}$|2\n",
    "$\\vdots$|$\\vdots$\n",
    "$X_{1}^{(k)}$|$k$\n",
    "$X_{2}^{(k)}$|$k$\n",
    "$\\vdots$|$\\vdots$\n",
    "$X_{n_k}^{(k)}$|$k$\n",
    "\n",
    "Presentamos a continuación la serie de pasos para realizar el ANOVA de un factor. Para ello, iremos ilustrando con el siguiente ejemplo:\n",
    "\n",
    "***\n",
    "***\n",
    "**Ejemplo 1.**\n",
    "\n",
    "Se realizó un estudio para investigar el efecto del CO2 sobre la tasa de crecimiento de Pseudomonas fragi (un corruptor de alimentos). Se cree que el crecimiento se ve afectado por la cantidad de CO2 en el aire.\n",
    "\n",
    "Para contrastarlo, en un experimento se administró CO2 a 5 presiones atmosféricas diferentes a 10 cultivos diferentes por cada nivel, y se anotó el cambio (en %) de la masa celular al cabo de una hora:\n",
    "\n",
    "Presión de CO2 (en atmósferas)|0.0|0.083|0.29|0.50|0.86\n",
    ":--|:--:|:--:|:--:|:--:|:--:\n",
    "|62.6|\t50.9|\t45.5\t|29.5|\t24.9\n",
    "|59.6|\t44.3|\t41.1\t|22.8|\t17.2\n",
    "|64.5|\t47.5|\t29.8\t|19.2|\t7.8\n",
    "|59.3|\t49.5|\t38.3\t|20.6|\t10.5\n",
    "|58.6|\t48.5|\t40.2\t|29.2|\t17.8\n",
    "|64.6|\t50.4|\t38.5\t|24.1|\t22.1\n",
    "|50.9|\t35.2|\t30.2\t|22.6|\t22.6\n",
    "|56.2|\t49.9|\t27.0\t|32.7|\t16.8\n",
    "|52.3|\t42.6|\t40.0\t|24.4|\t15.9\n",
    "|62.8|\t41.6|\t33.9\t|29.6|\t8.8\n",
    "\n",
    "La tabla la encuentras con el nombre [presiones.csv](https://github.com/scidatmath2020/Inferencia-Estad-stica-2022/blob/main/data/presiones.csv)\n",
    "\n",
    "***\n",
    "***\n",
    "\n",
    "## Pasos del ANOVA\n",
    "\n",
    "Recuerda tener tu tabla en formato largo. **A partir del Paso 2 es fundamental**.\n",
    "\n",
    "**Paso 1) Prueba de normalidad**\n",
    "\n",
    "$p$-valor|Decisión|Significado\n",
    ":--|:--|:--\n",
    "Pequeño|Rechazar $H_0$|Hay buena probabilidad de que **NO es gaussiana**\n",
    "Grande|Rechazar $H_1$|Hay buena probabilidad de que **SÍ es gaussiana**\n",
    "\n",
    "**Recuerda que para utilizarlos en R requieres la paquetería ``fBasics``**\n",
    "\n",
    "\n",
    "Prueba|Abreviatura|Tamaño de<br> muestra|Desventaja|Estadístico|**R**\n",
    ":--|:--:|:--:|:--|:--:|:--\n",
    "Kolmogorov-Smirnov-Lilliefors|KSL|$5\\le n$|No detecta bien las colas|D|``lillieTest()``\n",
    "Anderson-Darling|AD|$7\\le n\\le1000$|No trabaja bien en muestras grandes|A|``adTest()``\n",
    "Shapiro-Wilks|SW|En **R** admite $\\le5000$|Sensible a repeticiones|W|``shapiroTest()``\n",
    "D'Agostino|D|$20\\le n$||Omnibus|``dagoTest()``\n",
    "\n",
    "El $p-$ valor se obtiene con ``fBasics::nombreTest(muestra)@test$p.value``\n",
    "\n",
    "**Paso 2) Pruebas de homocedasticidad**\n",
    "\n",
    "El test mas usual es el de Bartlett. En **R** viene dado por ``bartlett.test`` (no necesita paqueterías) y su notación es \n",
    "\n",
    "``bartle.test(X~F,data=tabla)`` \n",
    "\n",
    "donde \n",
    "\n",
    "``X`` es el nombre de la columna de la característica $X$\n",
    "\n",
    "``F`` es el nombre de la columna factor\n",
    "\n",
    "``tabla`` es el nombre del dataframe\n",
    "\n",
    "La decisión se toma de la siguiente manera:\n",
    "\n",
    "$p$-valor|significado\n",
    ":--|:--\n",
    "Pequeño|Hay buena probabilidad de que las varianzas no sean iguales\n",
    "Grande|Hay buena probabilidad de que las varianzas sean iguales\n",
    "\n",
    "El $p$-valor se obtiene con ``bartlett.test(X~F,data=tabla)$p.value``\n",
    "\n",
    "**Paso 3) Realización del ANOVA**\n",
    "\n",
    "En **R**, la función para aplicar el ANOVA es ``aov()`` que no requiere instalar ninguna paquetería. Su sintaxis es \n",
    "\n",
    "``aov(X~F,data=tabla)`` \n",
    "\n",
    "donde \n",
    "\n",
    "``X`` es el nombre de la columna de la característica $X$\n",
    "\n",
    "``F`` es el nombre de la columna factor\n",
    "\n",
    "``tabla`` es el nombre del dataframe\n",
    "\n",
    "El $p$-valor se obtiene con ``summary(resultado_anova)[[1]]$\"Pr(>F)\"`` donde ``resultado_anova=aov(X~F,data=tabla)``\n",
    "\n",
    "La decisión se toma de la siguiente manera:\n",
    "\n",
    "$p$-valor|significado\n",
    ":--|:--\n",
    "Pequeño|Hay buena probabilidad de que las medias no sean iguales\n",
    "Grande|Hay buena probabilidad de que las medias sean iguales\n",
    "\n",
    "** Paso 4) Hallar las medias diferentes**\n",
    "\n",
    "Recuerda tener tu tabla en formato largo. \n",
    "\n",
    "Se consideran dos casos: si las subpoblaciones están o no balanceadas (es decir, el tamaño de cada submuestra es aproximadamente igual o no).\n",
    "\n",
    "En el caso de balanceados, utilizamos la prueba de Duncan. En el caso de no balanceadas, utilizamos la prueba de Tukey. Ambas están en la paquetería ``agricolae``. \n",
    "\n",
    "Muestra|Agrupamiento por grupos|Función|Grupos|Descripción\n",
    ":--|:--|:--|:--|:--\n",
    "Desbalanceada|Tukey|``HSD.test``|``HSD.test(resultado_anova,trt = columna de factor,group = FALSE)$comparison``|Devuelve los $p$-valores\n",
    "|||``HSD.test(resultado_anova,trt = columna de factor,group = TRUE)$groups``|Devuelve los grupos que tienen medias iguales\n",
    "Balanceada|Duncan|``duncan.test``|``duncan.test(resultado_anova,trt = columna de factor,group = FALSE)$comparison``|Devuelve los $p$-valores\n",
    "|||``duncan.test(resultado_anova,trt = columna de factor,group = TRUE)$groups``|Devuelve los grupos que tienen medias iguales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Práctica.\n",
    "\n",
    "La tabla [sismos.csv](https://drive.google.com/file/d/1QT-L1RLHlzO-vOJQFaI2lzdtC3vAKrnE/view?usp=sharing) consiste de la información del Servicio Sismológico Nacional de todos los sismos registrados desde el 1 de enero de 1900 hasta las 13:40pm (hora de CDMX) del 22 de septiembre de 2022. Realiza con ella un ANOVA para:\n",
    "\n",
    "1. Verificar que la magnitud promedio mensual de los sismos en México varía durante todo el año.\n",
    "\n",
    "2. Verificar que la magnitud promedio mensual de los sismos fuertes (arriba de 6.5) en México es la misma durante todo el año."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}