Planteamiento

• El objetivo de las pruebas de bondad de ajuste es contrastar las hipótesis \[ H_0: F = F_0 \qquad \text{vs.} \qquad F \neq F_0 \] donde \(F_0\) es una función de distribución completa o parcialmente conocida, a partir de una muestra aleatoria \(X_1, \ldots, X_n\) de \(F(x)\).

• La prueba de Pearson se utiliza directamente para distribuciones discretas con un número finito de saltos y se puede adaptar para distribuciones continuas o discretas con infinitos saltos (por ejemplo, las distribuciones geométrica o la de Poisson).

• El objetivo es comparar las frecuencias observadas en cada categoría contra las frecuencias esperadas bajo el supuesto que \(F = F_0\). Cuando la distribución no es discreta con un número finito de saltos se induce una discretización para formar categorías.

Procedimiento

• Si \(F_0\) es continua o tiene infinitos saltos, primero se determina el número de categorías a utilizar y los extremos de los intevalos que las definen.

• Si se eligen \(c\) categorías y constantes \[ -\infty = a_0 < a_1 < a_2 < \ldots < a_{c-1} < a_c = \infty \] de manera que se formen \(c\) intervalos donde el k-ésimo está dado por

\[ S_k = \left\lbrace \begin{array}{ll} (a_{k-1}, a_k], & \text{para } k = 1, \ldots, c-1. \\ (a_{c-1}, \infty), & \text{si } k = c. \end{array} \right. \]

• Una vez definidas las categorías a partir de los intervalos, se calculan las fecuencias observadas como \[ O_k = \sum_{i = 1}^n 1\{ X_i \in S_k \}, \qquad k = 1, \ldots, c. \]

• Las frecuencias esperadas se calculan como \[ E_k = nPr\{X \in S_k \,\vert\, H_0\} = n\left(F_0(a_k) - F_0(a_{k-1}) \right), \qquad k = 1, \ldots, c. \]

Procedimiento

• Si \(F_0\) es discreta con un número finito de saltos se puede pensar en una variable aleatoria \(X\) con un número finito de valores, por ejemplo, \(\{x_1, \ldots, x_c \}\). En este caso las frecuencias se calcular directamente como \[ O_k = \sum_{i=1}^n I\{ X_i = x_k \}, \qquad k = 1, \ldots, c. \]

• En este caso, \(F_0\) está caracterizada por los valores \(P\{X = x_k \} = p_{0k}\), por lo que las frecuencias esperadas se calculan directamente como \[ E_k = nP\{X = x_k \} = np_{0k}, \qquad k = 1, \ldots, c. \]

Procedimiento

• Para medir la discrepancia entre lo observado y lo esperado Pearson propuso utilizar el estadístico \[ X^2 = \sum_{k=1}^c \frac{(O_k - E_k)^2}{E_k} \]

• La distribución de \(X^2\) bajo \(H_0\) es aproximadamente \(\chi^2_{(c-1)}\), por lo que para rechazar \(H_0\) se compara \(X^2\) con los valores críticos esta distribución (véase Cramér, 1946: 416-419).

• La calidad de la aproximación depende de las frecuencias observadas. La recomendación es elegir el número de categorías y los extremos de los intervalos tales que \(O_e \geq 5\), \(k = 1, \ldots, c\).

• Otra alternativa es aproximar la distribución de \(X^2\) con simulación, ya que \(F_0\) está completamente específicada.

Aproximación de la distrabución de \(X^2\)

El vector de frecuencias observadas \(\mathbf{O} = (O_1, \ldots, O_c)\) sigue una distribución multinomial con parámetros \(n\) y vector de probabilidades \[ \mathbf{P}_0 = (Pr(X \in S_1\,\vert\, H_0), \ldots, Pr(X \in S_c\,\vert\, H_0)) \] o bien \[ \mathbf{P}_0 = (Pr\{X = x_1\}, \ldots, Pr\{X = x_c \}). \] Se puede utilizar el siguiente procedimiento para simular observaciones de la distribución verdadera de \(X^2\).

• Se genera un vector \(\mathbf{O}_1\) de la distribución \(Multinomial(n, \mathbf{P}_0)\).

• A partir de las frecuencias \(O_{1,1}, \ldots, O_{c,1}\) se calcula el estadístico \(X^2_1\).

• Lo anterior se repite un número grande de veces \(m\) para obtener \(X^2_1, \ldots, X^2_m\).

• La distribucion de frecuencias de los valores \(X^2_1, \ldots, X^2_m\) se puede utilizar para aproximar la distribución verdadera de \(X^2\), en cada caso en particular.

Ejemplo: distribución discreta

Una planta de jardín presenta dos variedades: una de flores rojas y hojas alargadas y otra de flores blancas y hojas pequeñas. El carácter color de las flores sigue una herencia intermedia y el carácter tamaño de la hoja presenta una dominancia del carácter alargado. Las proporciones fenotípicas esperadas en la segunda generación son las siguientes

En un experimento de cruza de especímenes de las dos varieades se obtuvieron 654 plantas en la segunda generación y se observaron los siguientes resultados.

Utilizar la prueba Ji-cuadrada para comprobar si los resultados del experimento son consistentes con la teoría.

• Como ya se tienen las frecuencias observadas, sólo resta calcular las frecuencias esperadas, éstas se calculan multiplicando el tamaño de muestra por la proporción esperada de cada fenotipo. En la siguiente tabla se muestran los resultados y las cantidades intermedias para calcular el estañistico \(X^2\).

donde E_k representa las frecuencias esperadas (\(E_k\)), Diff representa diferencia entre la frecuencia observada y la esperada (\(O_k-E_k\)), y Diff2_Ek representa \((O_k - E_k)^2/E_k\).

• De los resultados anteriores se sigue que el estadístico \(X^2\) vale

## [1] 2.768603

• Si se toma \(\alpha = 0.01\), el valor crítico para concluir que los resultados del experimento respaldan el modelo dado para las proporciones de los fenotipos, está dado por el cuantil \(0.99\) de la distribución \(\chi^2_{(5)}\), que vale 15.09. Los grados de libertad corresponden al número de categorías \(c = 6\) menos uno.

• Como \(X^2 = 2.77 < 15.1 = \chi^2_{(5)}(0.99)\) se concluye que el experimento valida el modelo dado para las proporciones de los fenotipos.

Ejemplo: distribución normal estándar

• Se utilizará la prueba \(\chi^2\) para determinar si una muestra con 47 obervaciones de una distribución \(t_{(1)}\), proviene de una distribución \(N(0,1)\).

set.seed(1010)
ejemplo1 <- rt(47, 1)

• Arbitrariamente se decide construir \(c = 4\) categorías que correspondan a los intervalos definidos por los cuartiles de la distribución \(N(0, 1)\), de manera que las fecuencias observadas sean iguales.

• Los extremos de los intervalos son:

a <- qnorm(c(0.25, 0.5, 0.75)); a

## [1] -0.6744898  0.0000000  0.6744898

• Se calculan las frecuencias observadas

O <- c()
O[1] <- sum(ejemplo1 <= a[1])
O[2] <- sum(ejemplo1 <= a[2]) - sum(ejemplo1 <= a[1])
O[3] <- sum(ejemplo1 <= a[3]) - sum(ejemplo1 <= a[2])
O[4] <- sum(ejemplo1 > a[3]); O

## [1] 13  7 11 16

• Las fecuencias esperadas son todas iguales, \(E_k = 0.25n = 11.75\).

• Se calcula el estadístico \(X^2\).

x2 <- sum((O-11.75)^2/11.75); x2

## [1] 3.638298

Solución aproximada asintótica

• Bajo \(H_0\) la distribución aproximada de \(X^2\) es \(\chi^2_{(3)}\), dado que se tomaron 4 categorias.

• Para decidir si se rechaza \(H_0\) podemos utilizar el p-value del estadístico calculado con los datos.

pchisq(x2, 3, lower.tail = F)

## [1] 0.3032633

• Del resultado anterior se concluye no rechazar \(H_0\), esto significa que no hay evidencia para rechazar que de que la verdadera distribución de los datos es \(N(0, 1)\).

Solución aproximada con simulación

• También podemos aproximar el p-value con simulación con el procedimiento descrito anteriormente.

• En R, se utiliza la función rmultinom para simular valores de la distribución multinomial.

• El siguiente código se utiliza para simular 75,000 observaciones de la distribución de \(X^2\) bajo \(H_0\).

o_sim <- rmultinom(75E3, 47, rep(0.25, 4))
x2_sim <- apply(o_sim, 2, function(x) sum((x-11.75)^2/11.75)) 

• Podemos comparar la distribución de los valores simulados contra la distribución \(\chi^2_{(3)}\).

• Para decidir si se rechaza \(H_0\) o no, se calcula el p-value a partir de los valores simulados.

mean(x2 < x2_sim)

## [1] 0.3018667

• Por lo anterior se concluye que no se rechaza \(H_0\).

Ejemplo: distribución exponencial

• Se utilizará la prueba \(\chi^2\) para determinar si una muestra con 23 obervaciones de una distribución \(Ga(3, 2)\), proviene de una distribución \(Exp(2/3)\).

ejemplo2 <- rgamma(23, 3, 2)

• Ahora elegimos forman categorías \(c = 5\) que correspondan a los intervalos definidos por los cuantiles 0.2, 0.4, 0.6 y 0.8 de la distribución \(Exp(2/3)\).

• Los extremos de los intervalos son:

a <- qexp(c(0.2, 0.4, 0.6, 0.8), 2/3); a

## [1] 0.3347153 0.7662384 1.3744361 2.4141569

• Se calculan las frecuencias observadas. Se puede ver que algunas de las fecuencias son menores a 5.

O <- c()
O[1] <- sum(ejemplo2 <= a[1])
O[2] <- sum(ejemplo2 <= a[2]) - sum(ejemplo2 <= a[1])
O[3] <- sum(ejemplo2 <= a[3]) - sum(ejemplo2 <= a[2])
O[4] <- sum(ejemplo2 <= a[4]) - sum(ejemplo2 <= a[3])
O[5] <- sum(ejemplo2 > a[4]); O

## [1] 0 4 8 9 2

• Las fecuencias esperadas son todas iguales, \(E_k = 0.20n = 4.6\).

• Se calcula el estadístico \(X^2\).

x2 <- sum((O-4.6)^2/4.6); x2

## [1] 12.86957

Solución aproximada asintótica

• Bajo \(H_0\) la distribución aproximada de \(X^2\) es \(\chi^2_{(4)}\), dado que se tomaron 5 categorias.

• Para decidir si se rechaza \(H_0\) podemos utilizar el p-value del estadístico calculado con los datos.

pchisq(x2, 4, lower.tail = F)

## [1] 0.01193102

• Del resultado anterior se concluye rechazar \(H_0\), esto significa la evidencia señala que la la verdadera distribución de los datos no es \(Exp(2/3)\).

Solución aproximada con simulación

• También podemos aproximar el p-value con simulación con el procedimiento descrito anteriormente.

• En R, se utiliza la función rmultinom para simular valores de la distribución multinomial.

• El siguiente código se utiliza para simular 75,000 observaciones de la distribución de \(X^2\) bajo \(H_0\).

o_sim <- rmultinom(75E3, 23, rep(0.2, 5))
x2_sim <- apply(o_sim, 2, function(x) sum((x-4.6)^2/4.6)) 

• Se puede comparar la distribución de los valores simulados contra la distribución \(\chi^2_{(4)}\). En este caso se nota como la aproximación asintótica no es tan buena como en el ejemplo anterior.

• Para decidir si se rechaza \(H_0\) o no, se calcula el p-value a partir de los valores simulados.

mean(x2 < x2_sim)

## [1] 0.00932

• De igual manera se concluye rechazar \(H_0\).

\(\chi^2\) para hipótesis compuestas

• Existe una versión, quizá más útil, de la prueba \(\chi^2\) cuando \(H_0\) es compuesta, es decir, cuando \(F_0\) tiene parámetros desconocidos.

• Suponer que se desea contrastar las hipótesis \(H_0: F(\cdot) = F_0(\cdot \,\vert\, \boldsymbol{\theta})\) vs. \(H_1: F \neq F_0(\cdot \,\vert\, \boldsymbol{\theta})\), donde \(\boldsymbol{\beta}\) es el vector de parámetros desconocidos.

• Al igual que cuando \(H_0\) es simple se elige un numéro de categorias \(c\) y de definen los intervalos \(S_1, \ldots, S_c\).

• La diferencia radica en que las probabilidades \(P_{0k}=Pr(X \in S_k \,\vert\, H_0)\) dependen del vector de parámetros desconocidos \(\boldsymbol{\theta}\) y a su vez, las fecuencias esperadas \(E_1, \ldots, E_c\) dependen de los parámetros desconocidos.

• La corrección consiste en estimar \(\boldsymbol{\theta}\) a partir de los datos y hacer un plug-in para estimar \(P_{0k}\) y \(E_k\).

• También existe una versión del Teorema de Pearson que justifica que la distribución asintótica del estadístico \(X^2\) cuando \(F_0\) tiene \(r\) parámetros desconocidos, es \(X^2_{(c-r-1)}\).

• Para estimar \(\boldsymbol{\theta}\) se recomienda utilizar los estimadores de mínima \(\chi^2\) que se definen precisamente como los estimadores que minimizan el valor del estadístico \(X^2\) calculado haciendo un plug-in para el cálculo de \(E_k\). El inconveniente es que estos estimadores no siempre son fáciles de calcular.

• En los libros de Conover (1999) y Villegas (2005) se pueden consultar más detalles sobre esta modificación de la prueba \(\chi^2\). En el libro de Crámer (1946) se puede consultar la demostración de ambos Teoremas de Pearson sobre las distribuciones asintóticas de los estadísticos \(\chi^2\).

Referencias

• Conover, W. J. (1999). Practical Nonparametric Statistics, 3ra edición. New York: John Wiley & Sons.

• Cramér, H. (1946). Mathematical Methods of Statistics. Princeton University Press.

• Villegas, M. Á. G. (2005). Inferencia estadística. Ediciones Díaz de Santos.