Preliminares
Las pruebas que se verán se basan en los rangos de las observaciones en lugar de las observaciones en sí. Esto traslada el problema de comparar las poblaciones a estudiar la distribución de la suma de particiones aleatorias del conjunto \(\{1, \ldots, n \}\).
Definición (Rangos sin empates) Sea un conjunto de datos \(x_1, \ldots, x_n\) sin empates, es decir, \(x_i \neq x_j\) si \(i \neq j\). Se define el rango \(x_i\), denotado por \(R(x_i)\), como la posición que ocupa esta observación cuando los datos se ordenan de menor a mayor. Expresado de otra forma, \(R(x_i) = j\) si y sólo si \(x_i = x_{(j)}\).
Ejemplo Si se tiene el conjunto de datos \(\{2.1, 8.3, 3.5, 4.2, 6.2\}\), entonces los rangos correspondientes son \(\{1, 5, 2, 3, 4 \}\). Lo anterior significa que \(2.1\) es la menor observación, \(3.5\) es la segunda menor observación, así sucesivamente hasta llegar a que \(8.3\) es la mayor de las observaciones.
En R se utiliza la función rank para obtener los rangos de un vector de datos.
x <- c(2.1, 8.3, 3.5, 4.2, 6.2)
rank(x)
## [1] 1 5 2 3 4
Definición (Rangos con empates) Sea un conjunto de datos \(x_1, \ldots, x_n\), tal que para algunos indices \(i \neq j\), \(x_i \neq x_j\). El rango de \(x_i\), denotado por \(R(x_i)\), se define a partir de la posición la posición que ocupa esta observación cuando los datos se ordenan de menor a mayor, cuando dos o más observaciones tienen el mismo valor que \(x_i\) (hay empate), \(R(x_i)\) se define como el promedio de las posiciones de las observaciones empatadas con \(x_i\).
La definición anterior se clarifica con un ejemplo.
Ejemplo Considerar el siguiente conjunto de datos \(\{4, 2, 1, 6, 2, 4, 5, 3, 4\}\). Se tienen dos empates, hay dos observaciones con el valor \(2\) y tres observaciones con el valor \(4\). Para calcular los rangos primero se asignan las posiciones de menor a mayor
##
## Datos 4 2 1 6 2 4 5 3 4
## Posicion 5 2 1 9 3 6 8 4 7
Ya que se tienen las posiciones, se rompen los empates asignando el promedio de las posiciones de las observaciones empatadas. Por ejemplo, los dos observaciones empatadas con el valor \(2\), tienen las posiciones \(2\) y \(3\), por lo que el rango que le corresponde a cada una es \(2.5\).
##
## Datos 4 2.0 1 6 2.0 4 5 3 4
## Rango 6 2.5 1 9 2.5 6 8 4 6
En R se utiliza función rank con la opción ties.method = 'average' para calcular los rangos con empates de acuercon con la definición dada anteriormente. La opción ties.method = 'first' devuelve las posciones a partir de las que se obtienen los rangos promediando.
Prueba de Wilcoxon
Planteamiento
Se tiene una muestra aleatoria bivariada \((X_1, Y_1), \ldots, (X_n, Y_n)\) de una población con distribución \(F_{X, Y}\). El objetivo es contrastar las hipótesis \[
H_0: F_X = F_Y \qquad \text{vs.} \qquad H_1: F_X \neq F_Y
\] donde \(F_x\) es la distibución verdadera de las \(X_i\) y \(F_Y\) es la distribución verdadera de las \(Y_j\). Si es razonable asumir que \(F_X\) y \(F_Y\) son (absolutamente) continuas, de manera que \(P\{X = Y\} = 0\), las hipótesis anteriores se pueden plantear como \[
H_0: P\{X > Y\} = 0.5 \qquad \text{vs.} \qquad H_1: P\{X > Y\} \neq 0.5
\]
Si adicionalmente se puede asumir que \(F_X\) y \(F_Y\) son iguales salvo por la localización, es decir, si existe una constante \(c\in\mathbb{R}\) tal que \(F_X(x) = F_Y(x + c)\), \(\forall x \in \mathbb{R}\), entonces las hipótesis anteriores se pueden plantear como \[
H_0: E(X) = E(Y) \qquad \text{vs.} \qquad H_1: E(X) \neq E(Y)
\] o bien, \[
H_0: Med(X) = Med(Y) \qquad \text{vs.} \qquad H_1: Med(X) \neq Med(Y)
\] donde \(Med(X)\) representa la mediana de \(X\).
Procedimiento
- Se calculan las diferncias \(D_i = Y_i - X_i\) y los rangos del sus valores absolutos \(\vert\,D_i \,\vert\), \(R_i\).
- El estadístico de la prueba de suma de rangos con signo de Wilcoxon es el siguiente \[
T = \sum_{i=1}^n \text{sgn}(D_i) R_i
\]
- El estadístico de la prueba de suma de rangos de Wilcoxon es el siguiente \[
V = \sum_{i=1}^n I\{D_i > 0 \} R_i
\]
- En ambos casos, la distribución de los estadísticos puede obtenerse fácilmente con simulación. Para la prueba de suma de rangos con signo, bajo \(H_0\) y en caso de que no haya empates, \(D_i \neq D_j\), entonces es posible simular valores de la distribución de \(T\) como sigue
- Para \(i = 1, \ldots, n\) se generan variables aleatorias \(W_i\) tales que \(P\{W_i = 1\} = P\{W_i = -1 \} = 0.5\).
- Se calcula \(t^{(1)} = \sum_{i=1}^n iW_i\).
- Lo anterior se repite un número grande de veces para obtener \(t^{(1)}, \ldots, t^{(m)}\).
- \(t^{(1)}, \ldots, t^{(m)}\) constituyen una realización de una muestra aleatoria de la distribución del estadístico \(T\).
- Para simular valores de la distribución de \(V\) bajo \(H_0\), en caso de que no haya empates, se procede como sigue:
- Para \(i = 1, \ldots, n\) se generan variables aleatorias \(Z_i\) tales que \(P\{Z_i = 1\} = P\{Z_i = 0\} = 0.5\).
- Se calcula \(v^{(1)} = \sum_{i=1}^n i Z_i\).
- Lo anterior se repite un número grande de veces para obtener \(v^{(1)}, \ldots, v^{(m)}\).
- \(v^{(1)}, \ldots, v^{(m)}\) constituyen una realización de una muestra aleatoria de la distribución de \(V\).
Ejemplo. Prueba de suma de rangos de Wilcoxon
En R se utiliza la función wilcox.test para ejecutar la prueba de suma de rangos de Wilcoxon. Para ilustrar su utilización se considerará en siguiente conjunto de datos simulados:
x <- rnorm(15, 5, 3)
y <- x + rnorm(15, 0, 1)
- Se realiza la prueba con la función
wilcox.test de R
wilcox.test(x, y, paired = TRUE, exact = TRUE)
##
## Wilcoxon signed rank test
##
## data: x and y
## V = 66, p-value = 0.7615
## alternative hypothesis: true location shift is not equal to 0
Los resultados indican que no hay evidencia para rechazar que
- Para hacer la prueba a mano primero se calcula el valor del estadístico \(V\) para este conjunto de datos.
d <- x - y
v_obs <- sum((d > 0) * rank(abs(d), ties.method = 'average'))
v_obs
## [1] 66
Se genera una muestra aleatoria de tamaño \(m = 50,000\) del estadístico \(V\).
v_sim <- c()
for (j in 1:5E4)
{
v_sim[j] <- sum(rbinom(15, 1, 0.5) * 1:15)
}
Con la muestra anterior se calcula el p-value de la prueba
2 * min(mean(v_obs <= v_sim), mean(v_obs >= v_sim))
## [1] 0.75712
- También se puede hacer a mano la prueba de suma de rangos con signo de Wilcoxon. Primero se calcula el estadístico \(T\) para este conjunto de datos.
t_obs <- sum(((d > 0) - (d < 0)) * rank(abs(d), ties.method = 'average'))
t_obs
## [1] 12
Se genera una muestra aleatoria de tamaño \(m = 50,000\) del estadístico \(T\).
t_sim <- c()
for (j in 1:5E4)
{
t_sim[j] <- sum(sign(runif(15, -1, 1)) * 1:15)
}
Con la muestra anterior se calcula el p-value de la prueba
2 * min(mean(t_obs <= t_sim), mean(t_obs >= t_sim))
## [1] 0.75696
- Se puede observar que el p-value de las dos pruebas es similar.
Prueba \(U\) de Mann-Whitney
Planteamiento
Se tienen dos muestras aleatorias independientes \(X_1, \ldots, X_n \sim F_1\) y \(Y_1, \ldots, Y_m \sim F_2\). El objetivo es contrastar las hipótesis. \[
H_0: F_1 = F_2 \qquad \text{vs.} \qquad H_1: F_1 \neq F_2
\]
Al igual que para la prueba de Wilcoxon, se puede consisdera los pares de hipótesis \[
H_0: P\{X > Y\} = 0.5 \qquad \text{vs.} \qquad H_1: P\{X > Y\} = 0.5
\] si ambas distribuciones son absolutamente continuas, o bien \[
H_0: E(X) = E(Y) \qquad \text{vs.} \qquad E(X) \neq E(Y)
\] en el caso en que las distribuciones sean idénticas salvo por su localización.
Procedimiento
Se calculan los rangos de la de las observaciones en la muestra combinada \(X_1, \ldots, X_n, Y_1, \ldots, Y_m\), denotados por \(R(X_i)\) o \(R(Y_j)\), según sea el caso. Estos rangos van de 1 hasta \(n+m\).
Se calculan las sumas de los rangos de las observaciones de cada muestra, \[
T_1 = \sum_{i=1}^nR(X_i) \qquad \text{y} \qquad T_2 = \sum_{j=1}^mR(Y_j)
\]
Una alternativa consiste en utilizar directamente como estadísticos a \(T_1\) o \(T_2\). Otra posibilidad consiste en calcular los estadísticos \[
U_1 = nm + \frac{n(n+1)}{2} - T_1 \qquad \text{y} \qquad U_2 = nm + \frac{m(m+1)}{2} - T_2
\] El estadístico \(U\) de Mann-Whitney se toma como el mínimo de las sumas anteriores. \[
U = \min\left( U_1, U_2 \right)
\]
La prueba consiste en rechazar \(H_0: F_1 = F_2\) para valores extremos de \(U\).
La distribución nula del estadístico \(U\).
Bajo \(H_0\), la distribución de las observaciones en cada muestra es la misma, por lo que los rangos de las observaciones se reparten aleatoriaente en cada muestra.
Si no hay empates, la distribución de \(T_1\) corresponde a la suma de \(n\) números seleccionados al azar de entre los primeros \(n+m\) enteros positivos. Por simetría, \(T_2\) se distribuye como la suma de \(m\) enteros seleccionados al azar de entre los primero \(n+m\) enteros positivos.
A partir de las distribuciones de \(T_1\) y \(T_2\) se puede obtener una aproximación de la distribución de \(U\), por ejemplo con simulación. Obviamente la distribución de \(U\) depende de \(n\) y \(m\).
Con el siguiente código en R se puede simular de la distribución de \(U\) con \(n=7\) y \(m = 11\).
n <- 7; m <- 11; N <- n+m
u_sim <- c()
for (i in 1:5E4)
{
t1 <- sum(sample(N, n))
t2 <- 0.5*N*(N+1) - t1
u1 <- n*m + 0.5*n*(n+1) - t1
u2 <- n*m + 0.5*m*(m+1) - t2
u_sim[i] <- min(u1, u2)
}
- La muestra anterior se puede utilizar para calcular cuantiles de la distribución de \(U\) o bien para calcular p-values.
quantile(u_sim, probs = c(0.01, 0.05, 0.10, 0.20))
## 1% 5% 10% 20%
## 11 17 20 24
Ejemplo. Prueba \(U\) de Mann-Whitney
En este ejemplo se comparan dos muestras de poblaciones similares, una \(N(5, 1)\) y otra \(N(5.5, 1)\). Para poder utilizar los valores de \(U\) generados previamente, se eligen \(n = 7\) y \(m = 11\).
- Primero se generan los datos.
x <- rnorm(7, 5, 1)
y <- rnorm(11, 5.5, 1)
- Después se calcula el estadístico \(U\) para este conjunto de datos.
rangos <- rank(c(x, y), ties.method = 'average')
t1 <- sum(rangos[1:7])
t2 <- sum(rangos[8:18])
u1 <- 7*11 + 0.5*7*8 - t1
u2 <- 7*11 + 0.5*11*12 - t2
u_obs <- min(u1, u2); u_obs
## [1] 28
- Finalmente se calcula el p-value de la prueba.
## [1] 0.37816
Con este resultado la conclusión es no rechazar \(H_0\), es decir, no hay evidencia suficiente para rechazar que las poblaciones de las que provienen las muestras son iguales.
- También es posible hacer la prueba con la función
wilcox.test de R, en este caso se debe especificar que las muestras no están relacionadas (o apareadas) con el parámetro paired = FALSE.
wilcox.test(x, y, paired = FALSE, exact = TRUE)
##
## Wilcoxon rank sum test
##
## data: x and y
## W = 28, p-value = 0.3749
## alternative hypothesis: true location shift is not equal to 0
