knitr::opts_chunk$set(message = FALSE, warning = FALSE)
library(kableExtra)
library(knitr)
library(DescTools)
library(BoutrosLab.plotting.general)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## Loading required package: cluster
## Loading required package: hexbin
## Loading required package: grid
## 
## Attaching package: 'BoutrosLab.plotting.general'
## The following object is masked from 'package:stats':
## 
##     dist
library(vcd)
## 
## Attaching package: 'vcd'
## The following object is masked from 'package:latticeExtra':
## 
##     rootogram

Pruebas de Bondad de Ajuste

Las pruebas de bondad de ajuste se aplican a datos de frecuencias en escalas nominales. Estos métodos nos permiten evaluar qué tan bien una muestra de observaciones de frecuencias de una población, se ajusta a una distribución hipotética de las frecuencias.

Un uso regular en biología de esta prueba, es la evaluación del resultado de cruces y las frecuencias esperadas de los fenotipos resultantes.

Ejemplo 1: Un geneticista de plantas evalúa el resultado de cruces en una población de plantas con flores amarillas y flores verdes: de un total de 100 semillas que sembró, obtuvo 84 plantas con flores amarillas y 16 con flores verdes. ¿Se ajusta este resultado a un sistema de alelos co-dominantes o a un sistema dominante-recesivo?

El estadístico chi-cuadrado (\(\chi^2\))

Este estadístico nos permite medir cuanto una muestra de una distribución observada, en escala nominal, difiere de una distribución hipotética (esperada).
\[\chi^2 = \sum_{i=1}^k\frac{(f_i - \hat f_i)^2}{\hat f_i}\] donde:

\(f_i\) es la frecuencia observada en la categoría i
\(\hat f_i\) es la frecuencia esperada en la categoría i
k es el número de categorías

Cálculos con el Ejemplo 1 A partir del Ejemplo 1 podemos establecer las siguientes hipótesis, asumiendo que esperamos una distribución según un sistema dominante-recesivo:
\(H_0\): la muestra proviene de una población con una proporción 3:1 de flores amarillas a flores verdes.
\(H_A\): la muestra proviene de una población que no tiene la proporción 3:1.
Tabla de frecuencias

amarillas verdes n
\(f_i\) 84 16 100
\((\hat f_i)\) (75) (25) 
#calculo de chi-cuadrado
k = 2
fm <- c(84,16)
fe <- c(75,25)
chic1 <- sum(((fm - fe)^2)/fe)
#valor critico de chi^2
nu = k - 1
chit <- qchisq(0.05,nu, lower.tail = FALSE)
chi <- data.frame(chic1,chit)
kable(chi, format = "markdown", col.names = c("chi^2 calculado", "chi^2 crítico"))
chi^2 calculado chi^2 crítico
4.32 3.841459

El valor mayor de \(\chi^2\) calculado con respecto al valor crítico, nos permite rechazar la hipótesis nula, y decir que nuestra muestra no se ajusta a la proporción 3:1.

Corrección para \(\chi^2\)

Los valores calculados de \(\chi^2\) provienen de datos discretos, y producirán valores discontinuos, sin embargo la distribución \(\chi^2\) es una función continua. Para evitar este error en el cálculo de \(\chi^2\), se utiliza una correción (Yates): \[\chi^2 = \sum_{i=1}^k\frac{(|f_i - \hat f_i| - 0.5)^2}{\hat f_i}\] Vamos a reevaluar nuestro ejemplo:

#calculo de chi-cuadrado
k = 2
fm <- c(84,16)
fe <- c(75,25)
chic1 <- sum(((abs(fm - fe) - 0.5)^2)/fe)
#chi2 de la distribución
nu = k - 1
chit <- qchisq(0.05,nu, lower.tail = FALSE)
#tabla resultados
chi <- data.frame(chic1,chit)
kable(chi, format = "markdown", col.names = c("chi^2 calculado", "chi^2 crítico"))
chi^2 calculado chi^2 crítico
3.853333 3.841459

Seguimos rechazando la \(H_0\), pero con un menor margen de diferencia. La corrección de Yates debe usarse cuando \(\nu\) = 1.

Prueba de bondad de ajuste (k = 2) con \(\chi^2\) en R

chisq.test(fm, p=fe, rescale.p = TRUE, correct = TRUE)
## 
##  Chi-squared test for given probabilities
## 
## data:  fm
## X-squared = 4.32, df = 1, p-value = 0.03767

El resultado es igual - rechazamos \(H_0\) - al tener un valor de p (error tipo I) menor de 0.05. Observa que aunque se le indica a la función que realice la corrección de Yates (correct = TRUE), no lo hace porque no es necesaria en este caso, como vimos anteriormente.

Prueba de bondad de ajuste con más de dos categorías

En el Ejemplo 2 vamos a probar la hipótesis genética de una distribución 9:3:3:1, para un sistema dialélico con dominancia. Flor amarilla (A) y semilla lisa (L) son los alelos dominantes, mientras que flor verde (V) y semilla rugosa (R) son los alelos recesivos.

#frec observada
fo <- c(152,39,53,6)
n = sum(fo)
#frec esperada
proporcion <- c(9,3,3,1)
np  = sum(proporcion)
fe <- proporcion*n/np
#crear tabla
fit <- t(fo)
fet <- t(fe)
frec <- rbind(fit,fet)
row <- c("fi","fe")
tabla <- cbind(row,frec)
knitr::kable(tabla, format = "markdown", col.names = c(" ","A-L","A-R","V-L","V-R"))
A-L A-R V-L V-R
fi 152 39 53 6
fe 140.625 46.875 46.875 15.625

Ahora calcularemos el valor de \(\chi^2\) y obtendremos el valor de \(\chi^2_{0.05,\nu}\) (usando qchisq o una tabla).

chick4 <- sum(((fo - fe)^2)/fe)
#chi2 de la distribución
k = 4
nu = k - 1
chit <- qchisq(0.05,nu, lower.tail = FALSE)
#tabla resultados
chi <- data.frame(chick4,chit)
kable(chi, format = "markdown", col.names = c("chi^2 calculado", "chi^2 crítico"))
chi^2 calculado chi^2 crítico
8.972444 7.814728

Podemos concluir que podemos rechazar la hipótesis nula, o sea que es poco probable que se cumpla la proporción 9:3:3:1.

Prueba de bondad de ajuste (k > 2) con \(\chi^2\) en R

chisq.test(fo, p=fe, rescale.p = TRUE)
## 
##  Chi-squared test for given probabilities
## 
## data:  fo
## X-squared = 8.9724, df = 3, p-value = 0.02966

Obtenemos el mismo resultado que con el procedimiento manual, rechazamos la \(H_0\).

El estadístico G (log-likelihood)

Cuando la cantidad de observaciones totales (n) no es muy grande, o cuando algún valor de \(|f_i - \hat f_i| \geq \hat f_i\), se recomienda el uso del estadístico G en lugar de \(\chi^2\). Este estadístico se calcula de la siguiente manera: \[G = 2\sum f_i ln \frac{f_i}{\hat f_i}\] o su equivalente: \[G = 2[\sum{f_i ln f_i} - \sum{f_i ln \hat f_i}]\]

Vamos a usar los datos del Ejemplo 2 para la prueba de hipótesis usando G.

#frec observada
fo <- c(152,39,53,6)
n = sum(fo)
#frec esperada
proporcion <- c(9,3,3,1)
np  = sum(proporcion)
fe <- proporcion*n/np
#crear tabla
fit <- t(fo)
fet <- t(fe)
frec <- rbind(fit,fet)
row <- c("fi","fe")
tabla <- cbind(row,frec)
knitr::kable(tabla, format = "markdown", col.names = c(" ","A-L","A-R","V-L","V-R"))
A-L A-R V-L V-R
fi 152 39 53 6
fe 140.625 46.875 46.875 15.625 
#calculo de G
G <- 2*sum(fo*log(fo/fe))
sprintf("G calculada: %.3f", G)
## [1] "G calculada: 10.833"
k = 4
nu = k - 1
chit <- qchisq(0.05,nu, lower.tail = FALSE)
sprintf("chi-cuadrado(0.05,3): %.3f", chit)
## [1] "chi-cuadrado(0.05,3): 7.815"

Con el estadístico G también podemos concluir que rechazamos con seguridad la hipótesis nula 9:3:3:1.

Prueba de bondad de ajuste con el estadístico G usando R

#pasar proporciones teóricas a probabilidades
prob1 <- proporcion/sum(proporcion)
# o pasar frecuencias esperadas a probabilidades
prob2 <- fe/sum(fe)
GTest(x=fo, p=prob1, correct="none")
## 
##  Log likelihood ratio (G-test) goodness of fit test
## 
## data:  fo
## G = 10.833, X-squared df = 3, p-value = 0.01267

Obtenemos los mismos resultados que con el cálculo manual.

Prueba de bondad de ajuste Kolmogorov-Smirnov

En las pruebas anteriores se comparan las frecuencias de la muestra y esperada, en una categoría nominal. La prueba de Kolmogorov-Smirnov nos permite comparar frecuencias observadas y esperadas, en escalas de proporciones (acumuladas), ordinales o intervalos. Esta prueba también se usa para probar si las distribuciones de dos muestras de valores continuos, son iguales o no (provienen o no de la misma población).

Datos en una escala de proporciones continuas

Para esta prueba se necesita calcular la frecuencia relativa observada acumulada de la variable observada: \[rel\ F_i = \frac{F_i}{n}\] donde:

\(F_i\) es la frecuencia acumulada observada hasta un valor de la escala continua, \(X_i\)
\(n\) es el total de observaciones realizadas a lo largo de la escala continua

La frecuencia relativa esperada acumulada se calcula a partir de la distribución esperada a lo largo de una escala continua: \[rel\ \hat F_i = \frac{X_i}{X_{máximo}}\]

El Ejemplo 3 se refiere a la distribución de mariposas en un árbol (Figura 1).

Figura 1. Distribución de mariposas a lo largo de un árbol. \(X_i\) es la distancia (m) desde el suelo a la que se encuentran \(f_i\) mariposas (frecuencia). Nota: solo se muestran algunos de los datos.

Datos de altura (m) a la que se observan las mariposas, \(X_i\), y número (frecuencia) de mariposas (\(f_i\)) observadas a la altura correspondiente. Queremos probar la hipótesis nula de que las mariposas se distribuyen uniformemente a lo largo de la altura del árbol (hasta 25 m), versus la hipótesis alterna de que las mariposas no se distribuyen uniformemente.

## [1] "Xi:"
##  [1]  1.4  2.6  3.3  4.2  4.7  5.6  6.4  7.7  9.3 10.6 11.5 12.4 18.6 22.3
## [1] "fi:"
##  [1] 1 1 1 1 1 2 1 1 1 1 1 1 1 1

Table 1. Frecuencia relativa acumulada de mariposas (\(rel\ F_i\)), frecuencia relativa acumulada de altura (\(rel\ \hat F_i\)), estadístico para la prueba de Kolmogorov-Smirnov: \(D_i = |rel\ F_i - rel\ \hat F_i|\).

#secuencia
i <- c(seq(1:14))
#frecuencia acumulada de mariposas
fiac <- cumsum(fi)
#frecuencia relativa acumulada de mariposas
firelac <- fiac/15
#frecuencia relativa acumulada de la altura
xirelac <- xi/25
#estadístico Di
di <- abs(firelac - xirelac)
#tabla de resultados
tablaKS <- data.frame(i, xi, fi, fiac, firelac, xirelac, di)
kable(tablaKS, format = "markdown", col.names = c("i","Xi","fi","Fi","rel F_i","rel F^_i","Di"))
i Xi fi Fi rel F_i rel F^_i Di
1 1.4 1 1 0.0666667 0.056 0.0106667
2 2.6 1 2 0.1333333 0.104 0.0293333
3 3.3 1 3 0.2000000 0.132 0.0680000
4 4.2 1 4 0.2666667 0.168 0.0986667
5 4.7 1 5 0.3333333 0.188 0.1453333
6 5.6 2 7 0.4666667 0.224 0.2426667
7 6.4 1 8 0.5333333 0.256 0.2773333
8 7.7 1 9 0.6000000 0.308 0.2920000
9 9.3 1 10 0.6666667 0.372 0.2946667
10 10.6 1 11 0.7333333 0.424 0.3093333
11 11.5 1 12 0.8000000 0.460 0.3400000
12 12.4 1 13 0.8666667 0.496 0.3706667
13 18.6 1 14 0.9333333 0.744 0.1893333
14 22.3 1 15 1.0000000 0.892 0.1080000

Para la prueba de hipótesis calculamos el valor máximo de \(D_i\), para compararlo con el valor crítico del estadístico: \(D_{0.05(2),n}\). D se obtiene a partir de una tabla o mediante la función ks.test.critical.value(n, p, alternative = “one-sided” (or “two-sided”)).

#valor máximo de Di
dmax <- max(di)
#D de la distribución
n = 15
p = 0.95
dcrit <- ks.test.critical.value(n, p, alternative = "two-sided")
tablaKStest <- data.frame(dmax,n, p, dcrit)
kable(tablaKStest, format = "markdown", col.names = c("D máximo", "n", "p", "D crítico"))
D máximo n p D crítico
0.3706667 15 0.95 0.3375964

El valor de D máximo es mayor que el D crítico, por lo tanto rechazamos la \(H_0\), y podemos decir que las mariposas no tienen una distribución uniforme a lo largo de la altura del árbol.

Tablas de Contingencia

La tabla de contingencia es parte de un procedimiento para probar la independencia entre variables nominales. La forma más común para probar la independencia es con el uso del estadístico \(\chi^2\). La forma más común de calcularlo es: \[\chi^2 = \sum_{i=1}^r \sum_{j=1}^c \frac{(f_{ij} - \hat f_{ij})^2}{\hat f_{ij}}\] donde:

\(f_{ij}\) es la frecuencia observada de la fila i, columna j
\(\hat f_{ij}\) es la frecuencia esperada de la fila i, columan j
c es el número de columnas
r es el número de filas (‘rows’)

Para calcular las frecuencias esperadas es necesario primero calcular las frecuencias totales de las filas (\(R_i\)) y de las columnas (\(C_j\)); para esto usamos las fórmulas: \[R_i = \sum_{j=1}^c f_{ij}\] \[C_j = \sum_{i=1}^r f_{ij}\] A partir de estos valores, podemos calcular cada frecuencia esperada: \[\hat f_{ij} = \frac{(R_i)(C_j)}{n}\] donde n es l número total de observaciones.
Una vez calculado el valor de \(\chi^2\), debemos compararlo con el valor crítico de la tabla (o función de R) para un nivel de significancia determinado (usualmente 0.05) y \(\nu = (r - 1)(c - 1)\).

Ejemplo 4. Queremos probar la independencia entre la característica color de pelo y sexo, en 300 humanos, seleccionando al azar 100 hombres y 200 mujeres.

color de pelo
sexo NEGRO MARRON RUBIO ROJO total
HOMBRE 32 43 16 9 100 = \(R_1\)
MUJER 55 65 64 16 200 = \(R_2\)
total 87 = \(C_1\) 108 = \(C_2\) 80 = \(C_3\) 25 = \(C_4\) 300 = \(n\)

La forma más general de plantear la hipótesis de independencia es: \(H_0\): el color del pelo es independiente del sexo en la población muestreada, y \(H_A\): el color del pelo no es independiente del sexo en la población muestreada.

Cálculo manual de \(\chi^2\)

Calculamos las frecuencias esperadas para calcular luego el estadístico \(\chi^2\).

Ri <- c(100,200)
Ci <- c(87,108,80,25)
mfe <- (Ri%*%t(Ci))/300
fe <- as.data.frame(mfe)
sexo <- c("Hombre","Mujer")
colnames(fe) <- c("Negro","Marrón","Rubio","Rojo")
tablafe <- cbind(sexo,fe)
kable(tablafe, format = "markdown", caption = "Frecuencias Esperadas")
sexo Negro Marrón Rubio Rojo
Hombre 29 36 26.66667 8.333333
Mujer 58 72 53.33333 16.666667 
#matriz de frecuencia esperada es mfe
#definimos matriz frecuencias observadas
mfo <- matrix(c(32,43,16,9,55,65,64,16), ncol = 4, byrow = TRUE)
#cálculo de chi-cuadrado
chicu2 = 0
for(i in 1:2){
  for(j in 1:4){
    chicu1 = sum(((mfo[i,j] - mfe[i,j])^2)/mfe[i,j])
      chicu2 = chicu2 + chicu1
  }
}
#chi-cuadrado critico de función (o tabla)
r = 2
c = 4
nu = (r-1)*(c-1)
chicri <- qchisq(0.05,nu, lower.tail = FALSE)
chitest <- data.frame(chicu2,nu,chicri)
colnames(chitest) <- c("Chi^2 calculado", "grados de libertad", "Chi^2 crítico")
kable(chitest, format = "markdown")
Chi^2 calculado grados de libertad Chi^2 crítico
8.987184 3 7.814728

El valor calculado de \(\chi^2\) resultó mayor que el valor crítico, por lo tanto rechazamos la hipótesis nula de independencia, y podemos esperar que el color del pelo no sea independiente del sexo.

Prueba de \(\chi^2\) utilizando R

#data frame desde matriz de datos
fo <- as.data.frame(mfo)
chisq.test(fo)
## 
##  Pearson's Chi-squared test
## 
## data:  fo
## X-squared = 8.9872, df = 3, p-value = 0.02946

Gráfica de mosaico para tablas de contingencia

Ahora bien, la prueba estadística solo nos indica si se acepta o rechaza la hipótesis nula, pero no nos específica dónde están las principales discrepancias entre las frecuencias, por efecto de la no-independencia entre los variables nominales. Una gráfica de mosaico puede ayudarnos a identificar las categorías entre las que se producen las mayores discrepancias.

mosaic(mfo, shade = TRUE, legend = FALSE, labeling_args = list(set_varnames = c(A = "Sexo", B = "Color Pelo")), set_labels = list(A = c("Hombre","Mujer"), B = c("Negro","Marrón","Rubio","Rojo")))