\documentclass[draft]{article}
\usepackage{beamerarticle}
\usepackage[utf8]{inputenc}
\usepackage[spanish]{babel}
\usepackage{graphicx}
\newcommand{\R}{{\mathbb R}}
\theoremstyle{definition}
\newtheorem{definicion}{Definición}
\newtheorem{teorema}{Teorema}
\newtheorem{ejemplo}{Ejemplo}
\renewcommand{\(}{\left(}  \renewcommand{\)}{\right)}
\renewcommand{\le}{\leqslant}
\renewcommand{\ge}{\geqslant}
\newcommand{\sigue}{\,\hookrightarrow\,}
\newcommand{\sigueHnula}{\stackrel{H_0}{\sigue}}
\newcommand{\unif}[2]{\mathcal U\(#1,#2\)}
\title{Contrastes múltiples}
\begin{document}
\maketitle
\begin{definicion}
  Llamaremos contraste múltiple a un contraste simultáneamente
  sobre varios parámetros.
\end{definicion}
\begin{ejemplo}
  El análisis de varianza contrasta la hipótesis nula
  $H_0\colon \mu_1=\dots=\mu_I$ frente a
  $H_1\colon \exists\,i,j,\;\mu_i\neq\mu_j$, por lo que es un
  contraste múltiple.
\end{ejemplo}
\begin{ejemplo}
  El contraste de homoscedasticidad de Bartlett contrasta
  $H_0\colon \sigma_1=\dots=\sigma_I$ frente a $H_1\colon
  \exists\,i,j,\;\sigma_i\neq\sigma_j$, luego también es un contraste múltiple.
\end{ejemplo}

La hipótesis nula de tales contrastes puede formularse como intersección de
hipótesis nulas parciales.

\begin{ejemplo}
  La hipótesis nula del contraste de Bartlett es $H_0=\bigcap_{i,j}H_0^{ij}$
  donde $H_0^{ij}\colon\sigma_i=\sigma_j$.  Sería la intersección de
  $k=\binom{I}{2}$ hipótesis nulas.
\end{ejemplo}

La hipótesis nula global $H_0$ se rechazaría al rechazar alguna de las
hipótesis parciales $H_0^{i}$, luego la región crítica global sería
\begin{equation}
\text{RC} = \bigcup_{i=1}^k \text{RC}_{i} \label{regcrit}
\end{equation}
donde $RC_{i}$ es la región crítica asociada al $i$-ésimo contraste parcial.

Para garantizar que
$$\Pr[\text{RC}\mid H_0]\le\alpha$$
se pueden considerar varios criterios.
\section{Bonferroni}
Aplicando la desigualdad de Boole,
$$ \Pr\(\bigcup_{i=1}^k A_i\)\le\sum_{i=1}^k\Pr A_i $$
a las regiones críticas de \eqref{regcrit}
y resolviendo cada contraste parcial a nivel de significación
$\alpha \over k$ se tiene
$$
\Pr(\text{RC}\mid H_0) = \Pr\(\bigcup_{i=1}^k \text{RC}_i\;\middle|\;H_0\)
\le\sum_{i=1}^k\Pr(\text{RC}_i\mid H_0^i)\le\sum_{i=1}^k{\alpha\over k}=\alpha
$$
\begin{ejemplo}
  En un contraste para comparar la igualdad de tres medias,
  \begin{eqnarray*}
    H_0\colon&&\mu_1=\mu_2=\mu_3 \\
    H_1\colon&&\mu_1\neq\mu_2\vee\mu_1\neq\mu_3\vee\mu2\neq\mu_3
  \end{eqnarray*}
  para conseguir un nivel de significación global $\alpha=0.05$
  se realizaría cada contrastre $H_0^{ij}\colon \mu_i=\mu_j$,
  $H_1^{ij}\colon \mu_i\neq\mu_j$ a nivel de significación
  ${\alpha\over3}={0.05\over3}=0.01666...$
\end{ejemplo}

El método de Bonferroni es muy popular porque
\begin{itemize}
\item Es muy fácil de aplicar.
\item Se puede usar siempre.
\end{itemize}
Su inconveniente fundamental es que el tamaño del contraste global
puede estar muy alejado de $\alpha$, con la consiguiente pérdida de potencia.
\section{\v Sid\'ak}
Cuando las regiones críticas de los contrastes parciales se basan en
estadísticos independientes, entonces
\begin{multline*}
\alpha=  \Pr(\text{RC}\mid H_0) =
1-\Pr(\text{RA}\mid H_0) =
1-\Pr\(\bigcap_{i=1}^k \text{RA}_i\;\middle|\;H_0\)=\\
=1-\prod_{i=1}^k\Pr(\text{RA}_i\mid H_0^i)=1-\prod_{i=1}^k(1-\alpha_i)
\end{multline*}
luego bastaría tomar como nivel de significación de cada
contraste parcial
$$\alpha_i=1-\sqrt[k]{1-\alpha}$$

Aunque esta corrección es menos conservadora que la de Bonferroni
y, por tanto, más potente, suele serlo por poco.
\begin{ejemplo}
$$
\left.
  \begin{array}{l}
    \alpha=0.05\\k=10
  \end{array}
\right\}
\implies
\left\{
  \begin{array}{l}
    \text{Bonferroni } \alpha_i=0.005\\
    \text{\v Sid\'ak } \alpha_i=0.005\,116
  \end{array}
\right.
$$
\end{ejemplo}
\section{Holm}
Este criterio mejora la corrección de Bonferroni, pues
\begin{itemize}
\item También puede usarse siempre (no requiere independencia).
\item En caso de rechazarse alguna hipótesis nula parcial,
  puede rechazar más hipótesis que el de Bonferroni.
\end{itemize}

Considérese la siguiente notación:
\begin{description}
\item[$p_{i}$] p-valor del contraste $H_0^i/H_1^i$, con $i=1,\dots,k$.
\item[$p_{(i)}$] p-valor ordenado, tal que $p_{(1)}\le\dots\le p_{(k)}$.
\item[$H_0^{(i)}$] hipótesis nula del contraste parcial con p-valor $p_{(i)}$.
\end{description}

El criterio de Holm consiste en
rechazar $H_0^{(1)}$, \dots, $H_0^{(m)}$ 
y no rechazar $H_0^{(m+1)}$, \dots, $H_0^{(k)}$
si
$$
\forall\,i\le m\quad
p_{(i)}\le\frac{\alpha}{k-i+1}
\qquad\wedge\qquad
p_{(m+1)}>\frac{\alpha}{k-m}
$$

Vamos a comprobar que el procedimiento de Holm mantiene un nivel de
significación global $\alpha$:
\begin{itemize}
\item Sea $I_0$ = $\{i\mid H_0^{(i)}\text{ es cierta}\}$
\item Sea $k_0=\#I_0$ el número de hipótesis nulas ciertas.
\item Supongamos que $H_0^{(j)}$ es la primera hipótesis nula cierta
  pero rechazada.
\item $H_0^{(1)}$, \dots, $H_0^{(j-1)}$ son falsas.
\item $j-1\le k-k_0$
\item $k_0\le k-j+1$
\item $j\le m \implies p_{(j)}\le\frac{\alpha}{k-j+1}\le\frac{\alpha}{k_0}$
\item Por tanto, si rechazamos una hipótesis nula verdadera, tiene que
  haber una hipótesis nula verdadera con p-valor $\le\frac\alpha{k_0}$.
\item Por tanto,
\begin{eqnarray*}
\Pr[\text{error global de tipo I}]
  &=&\Pr\(\text{rechazar alguna }H_0^{(i)}\text{ verdadera}\) \\
  &\le&\Pr\(\exists\,i\in I_0, p_i\le\frac\alpha{k_0}\) \\
  &\le&\sum_{i\in I_0}\Pr\(p_i\leq\frac\alpha{k_0}\)\\
&=& \sum_{i\in I_0}\Pr\(\text{RC}_i\text{ a nivel }\frac\alpha{k_0}\)\\
  &\le& \sum_{i\in I_0}\frac\alpha{k_0}=\alpha
\end{eqnarray*}
porque $\text{p-valor}\le a$ si y sólo si la muestra pertenece a
la región crítica construida a nivel de significación $a$. De hecho,
si la población es continua, $p_i\sigueHnula\unif{0}{1}$.

\end{itemize}
\section{R}
A veces, en vez de dividir el $\alpha$ global, se expresan los
p-valores de los contrastes sencillos multiplicados por el correspondiente
coeficiente.  Por ejemplo, al usar en R las funciones
\texttt{pairwise.prop.test}, \texttt{pairwise.t.test} y
\texttt{pairwise.}\-\texttt{wilcox.}\-\texttt{test}:
\begin{verbatim}
> summary (cw <- chickwts)

     weight             feed   
 Min.   :108.0   casein   :12  
 1st Qu.:204.5   horsebean:10  
 Median :258.0   linseed  :12  
 Mean   :261.3   meatmeal :11  
 3rd Qu.:323.5   soybean  :14  
 Max.   :423.0   sunflower:12  


> pairwise.t.test (cw$weight, cw$feed, "none") # p-valores crudos

	Pairwise comparisons using t tests with pooled SD 

data:  cw$weight and cw$feed 

          casein  horsebean linseed meatmeal soybean
horsebean 2.1e-09 -         -       -        -      
linseed   1.5e-05 0.01522   -       -        -      
meatmeal  0.04557 7.5e-06   0.01348 -        -      
soybean   0.00067 0.00032   0.20414 0.17255  -      
sunflower 0.81249 8.2e-10   6.2e-06 0.02644  0.00030

P value adjustment method: none 


> pairwise.t.test (cw$weight, cw$feed, "bonferroni")

	Pairwise comparisons using t tests with pooled SD 

data:  cw$weight and cw$feed 

          casein  horsebean linseed meatmeal soybean
horsebean 3.1e-08 -         -       -        -      
linseed   0.00022 0.22833   -       -        -      
meatmeal  0.68350 0.00011   0.20218 -        -      
soybean   0.00998 0.00487   1.00000 1.00000  -      
sunflower 1.00000 1.2e-08   9.3e-05 0.39653  0.00447

P value adjustment method: bonferroni 


> pairwise.t.test (cw$weight, cw$feed) # "holm" por omisión

	Pairwise comparisons using t tests with pooled SD 

data:  cw$weight and cw$feed 

          casein  horsebean linseed meatmeal soybean
horsebean 2.9e-08 -         -       -        -      
linseed   0.00016 0.09435   -       -        -      
meatmeal  0.18227 9.0e-05   0.09435 -        -      
soybean   0.00532 0.00298   0.51766 0.51766  -      
sunflower 0.81249 1.2e-08   8.1e-05 0.13218  0.00298

P value adjustment method: holm 
\end{verbatim}

Se puede solicitar la corrección de p-valores directamente:
\begin{verbatim}
> p.adjust (c(0.02,0.03), "none")
[1] 0.02 0.03
> p.adjust (c(0.02,0.03), "bonferroni")
[1] 0.04 0.06
> p.adjust (c(0.02,0.03), "holm")
[1] 0.04 0.04

> p.adjust (c (0.001, 0.005, 0.01, 0.015, 0.02, 0.04), "bonf")
[1] 0.006 0.030 0.060 0.090 0.120 0.240
> p.adjust (c (0.001, 0.005, 0.01, 0.015, 0.02, 0.04), "holm")
[1] 0.006 0.025 0.040 0.045 0.045 0.045
\end{verbatim}
\end{document}
