COMPARACIÓN ESTADÍSTICA DE MODELOS PREDICTIVOS

INTRODUCCIÓN

     El uso correcto de técnicas de selección de algoritmos es vital en trabajos de investigación, así como en aplicaciones de negocio. Acertar en la predicción de datos futuros es uno de los principales desafíos a la hora de desarrollar nuevos modelos predictivos o algoritmos.

     Cuando disponemos de resultados de dos o más modelos de machine learning (aplicados sobre mismos o distintos conjuntos de datos) y necesitamos compararlos formalmente para determinar cuál se comporta mejor ante un determinado problema, no sería formalmente válido escoger aquel que mejor resultado o métrica haya dado. Hay que demostrar que las diferencias son estadísticamente significativas, y no debidas al azar o ruido en los datos, de lo contrario no podemos asumir que los modelos en cuestión se comportan de forma distinta.

Desde el punto de vista estadístico, se establece una hipótesis nula a refutar:

\(\small{H_0}\): los modelos presentan el mismo rendimiento, dicho de otra forma, las métricas escogidas para evaluar su desempeño son iguales.

     Para dicho fin, tenemos a nuestra disposición un conjunto de técnicas o test estadísticos que podemos aplicar a los resultados de la calibración de nuestros modelos. En este sentido, podemos hablar de test estadísticos paramétricos y no paramétricos:

• Test paramétricos: asumen que los datos se distribuyen acorde a una determinada distribución estadística conocida (normal, Bernoulli, etc.). Además, son sensibles a outliers y suposiciones de independencia, normalidad y homocedasticidad. Por tanto, los datos deben cumplir ciertas condiciones para que los resultados de estos test sean fiables.

• Test no paramétricos: no asumen ninguna distribución en los datos. Los test paramétricos suelen preferirse sobre los no paramétricos ya que suelen caracterizarse por un mayor poder estadístico, esto es, son capaces con mayor probabilidad de detectar una diferencia estadísticamente significativa si esta existe.

     En función del problema, podemos recurrir a diferentes técnicas de remuestreo (validación cruzada de \(\small{k}\) pliegues, LOOCV, etc.) para calibrar nuestros modelos con datos de entrenamiento, así como diferentes métricas para determinar la eficiencia de los mismos.

     Un aspecto importante a tener en cuenta a la hora de determinar el tipo de test estadístico apropiado a aplicar es si contamos con medidas pareadas, es decir, métricas de los distintos modelos sobre el mismo conjunto de datos, o medidas sobre muestras independientes. En función de esto, podemos encontrarnos con distintos escenarios:

• Dos modelos en un mismo dominio o conjunto de datos.

• Varios modelos en un mismo dominio o conjunto de datos.

• Dos modelos en dominios o conjuntos de datos distintos.

• Varios modelos en dominios o conjuntos de datos distintos.

Los test aplicables para cada caso se recogen en el siguiente esquema:

     Así pues, la elección de uno u otro tipo de test depende, no solo de las necesidades del investigador, sino de las características del conjunto de datos y del tipo de métricas que queramos comparar.

Nota: En este documento se muestran ejemplos para casos con muestras dependientes.

DOS MODELOS EN UN MISMO DOMINIO

     Supongamos que contamos con un modelo clasificador A y un modelo clasificador B, y que para la calibración se ha aplicado una validación cruzada de \(\small{k}\) pliegues. Por cada pliegue, se utiliza un subconjunto de entrenamiento y validación distinto, a partir de los cuales obtenemos \(\small{k}\) medidas o métricas, que serán las que compararemos estadísticamente. Los test de este apartado se ejemplificarán en base a este ejemplo.

Test t de Student pareado

     Permite determinar si la media de las diferencias entre dos muestras pareadas es significativa, es decir, diferente de 0.

HIPÓTESIS

\(\small{H_0}\): \(\small{µ_A}\) - \(\small{µ_B}\) = 0 (La diferencia entre las medias pareadas de ambos grupos es igual a 0).

\(\small{H_1}\): \(\small{µ_A}\) - \(\small{µ_B}\) ≠ 0 (La diferencia entre las medias pareadas de ambos grupos es distinta de 0).

CONDICIONES

• Mismas muestras para cada grupo.

• Aleatoriedad de la muestra: las muestras a partir de las cuales se calculan las medias deben haber sido seleccionadas de forma independiente, siendo estas representativas.

• Normalidad: Las métricas de ambos modelos deben seguir una distribución normal. Entre los test comúnmente aplicables para comprobar este requerimiento se encuentran: Kolmogorov-Smirnov (recomendado con \(\small{n > 50}\)), Saphiro Wilk (recomendado con \(\small{n < 50}\)). Si esta condición no se cumple, el test \(\small{t}\) de Student puede ser robusto con conjuntos de test con más de 30 muestras. En el caso de aplicar validación cruzada con \(\small{k}\) pliegues, necesitaríamos un conjunto con \(\small{k}\) x 30 muestras. Otra opción sería aplicar esta validación cruzada con repetición para obtener al menos 30 medidas por pliegue.

ESTADÍSTICO

     Para el ejemplo comentado anteriormente, para cada observación \(\small{i}\) del conjunto de datos \(\small{n}\), el test calculará las \(\small{m}\) diferencias (correspondientes al número de pliegues de validación cruzada) entre el modelo A y el modelo B:

La ecuación para calcular el estadístico \(\small{t}\) sería:

donde

Este estadístico sigue una distribución \(\small{t}\) de Student con \(\small{m-1}\) grados de libertad. Con una significancia de \(\small{α/2}\) (normalmente 0,05), no se rechaza la hipótesis nula \(\small{H_0}\) si:

     Con este test podemos saber si la diferencia entre modelos es significativa, pero no nos aporta información sobre cuán importante es la diferencia, si es que existe. Para ello se calcula el estadístico d de Cohen, el cual es una medida del tamaño del efecto en base a la diferencia estandarizada de dos medias:

     El valor de este estadístico corresponde con las desviaciones típicas de diferencia entre los dos grupos. Una escala aproximada para interpretar el tamaño del efecto según este estadístico sería:

\(\small{d_{cohen}}\) ≈ 0,2-0,3 (pequeño)

\(\small{d_{cohen}}\) ≈ 0,5 (mediano)

\(\small{d_{cohen}}\) ≈ 0,8 (grande)

Consideraciones:

     Aunque este test es potente en base al uso de una validación cruzada con 10 pliegues, puede en este caso presentar un Error de Tipo I (detectar diferencias cuando no las hay) alto. Se recomienda en casos en los que el Error de Tipo II (fallar en detectar una diferencia real entre modelos) es más importante.

Test de McNemar

     El test de McNemar se presenta como alternativa no paramétrica al test \(\small{t}\) de Student para determinar la dependencia de datos categóricos pareados (clasificadores). Se aplica a tablas de contingencia 2x2 donde se incluye el número de muestras correcta e incorrectamente clasificadas por dos métodos o modelos considerados, para determinar si las frecuencias marginales de las filas y columnas de la tabla sin iguales (homogeneidad marginal):

	Modelo B (positivo)	Modelo B (negativo)	Total fila
Modelo A (positivo)	\(\small{n_{11}}\)	\(\small{n_{10}}\)	\(\small{(n_{11}+n_{10})}\)
Modelo A (negativo)	\(\small{n_{01}}\)	\(\small{n_{00}}\)	\(\small{(n_{01}+n_{00})}\)
Total columna	\(\small{(n_{01}+n_{01})}\)	\(\small{(n_{10}+n_{00})}\)	\(\small{n}\)

HIPÓTESIS

La hipótesis nula establece que ambos clasificadores tienen el mismo ratio de error:

\(\small{H_0}\): \(\small{n_{01}}\) = \(\small{n_{10}}\)

Nota: Para aplicar el test, comprobar que \(\small{n_{01}}\) + \(\small{n_{10}}\) > 20.

ESTADÍSTICO

     El estadístico, que se ajusta a una distribución Chi-cuadrado (\(\small{\chi^2}\)), se estima a partir de la ecuación con 1 grado de libertad

     El cálculo se centra en parejas discordantes (muestras clasificadas correctamente por un modelo, pero incorrectamente por el otro), siendo este el ratio de la diferencia al cuadrado de frecuencias discordantes en relación al total de frecuencias discordantes. Con una significancia \(\small{α}\), la hipótesis nula no se rechaza si \(\small{M ≤ {\chi^2}_{α,1}}\).

Ejemplo en R

Caso paramétrico

     Imaginemos que disponemos de un conjunto de datos sobre los que hemos entrenado dos clasificadores (C5.0 y Random Forest) mediante el método de validación cruzada de 10 pliegues y como métrica el índice Kappa. Por tanto, disponemos de 10 valores del índice Kappa para cada pliegue por clasificador:

metricas <- data.frame(C50 = c(0.9073259, 0.6822547, 0.6925111, 0.2719983, 
                               0.9073259, 0.0925111, 0.0925111, 0.8572170,
                               0.9106929, 0.0925111),
                       Rf = c(0.6541556, 0.7521948, 0.8056707, 0.2816065, 
                              0.9874889, -0.1236222, 0.3208222, 0.7521948, 
                              0.8056707, 0.1236222))

metricas

##          C50         Rf
## 1  0.9073259  0.6541556
## 2  0.6822547  0.7521948
## 3  0.6925111  0.8056707
## 4  0.2719983  0.2816065
## 5  0.9073259  0.9874889
## 6  0.0925111 -0.1236222
## 7  0.0925111  0.3208222
## 8  0.8572170  0.7521948
## 9  0.9106929  0.8056707
## 10 0.0925111  0.1236222

     Nos encontramos ante un caso de comparación de dos clasificadores ante un mismo dominio, ya que los dos se han calibrado sobre los mismos pliegues de datos (debemos asegurarnos de haber establecido una semilla para ello).

     Para comprobar si podemos aplicar un test \(\small{t}\) de Student de medidas pareadas, vamos a aplicar un test de Saphiro-Wilk para ver primero si se cumple la normalidad sobre las diferencias entre las métricas obtenidas en cada pliegue:

# Test Shapiro-Wilk
shapiro.test(metricas$C50 - metricas$Rf)

## 
##  Shapiro-Wilk normality test
## 
## data:  metricas$C50 - metricas$Rf
## W = 0.95884, p-value = 0.7725

     En base al resultado podemos afirmar, con un 95% de confianza, que los datos siguen una distribución normal (\(\small{W}\) = 0,958, p-value = 0,772).

     A continuación pasamos a comprobar el supuesto de homocedasticidad aplicando el test de Bartlett, que es más potente que el test de Levene bajo supuestos de normalidad:

# Metricas stack: concatemanos los valores de las metricas de cada modelo en un
# unico vector, con un segundo vector indicando de qué modelo procede el dato
metricas_stack <- stack(list(C50 = metricas$C50, Rf = metricas$Rf))
head(metricas_stack)

##      values ind
## 1 0.9073259 C50
## 2 0.6822547 C50
## 3 0.6925111 C50
## 4 0.2719983 C50
## 5 0.9073259 C50
## 6 0.0925111 C50

# Test de Bartlett
bartlett.test(values ~ ind, data = metricas_stack)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  values by ind
## Bartlett's K-squared = 0.0037017, df = 1, p-value = 0.9515

     El test confirma, con un 95% de confianza, que las metricas presentan homogeneidad de varianzas en cada modelo.

     Una vez cumplidas las condiciones para aplicar un test \(\small{t}\) de Student de medidas pareadas, pasamos a comprobar si las diferencias de rendimiento entre los dos modelos es significativa:

# Test t de Student para medidas pareadas
t.test(metricas$C50, metricas$Rf, paired = TRUE)

## 
##  Paired t-test
## 
## data:  metricas$C50 and metricas$Rf
## t = 0.30588, df = 9, p-value = 0.7667
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.09404931  0.12346029
## sample estimates:
## mean of the differences 
##              0.01470549

     No podemos rechazar, con un 95% de confianza, que la diferencia entre las medias de las métricas entre los modelos Rf y C50 son distintas a 0. Es decir, el rendimiento del modelo Rf no es significativamente distinto al C50 (\(\small{|T|}\) = 0,305, p-value = 0,766). De haber sido significativa, el signo del mean of the differences nos indicaría cual de los dos sería el mejor (Rf).

Caso no paramétrico

     Seguimos con el supuesto de contar con conjunto de métricas (índice Kappa) obtenidas de la calibración de dos modelos clasificadores (C5.0 y Random Forest) mediante el método de validación cruzada de 10 pliegues:

metricas <- data.frame(C50 = c(0.9333333, 0.8750000, 0.8125000, 0.8125000,
                               0.9333333, 0.8000000, 0.6000000, 0.9375000,
                               0.7500000, 0.8000000),
                       Rf = c(0.9333333, 0.8750000, 0.8750000, 0.8750000,
                              0.9333333, 0.8000000, 0.5333333, 0.9375000,
                              0.9375000, 0.7333333))

metricas

##          C50        Rf
## 1  0.9333333 0.9333333
## 2  0.8750000 0.8750000
## 3  0.8125000 0.8750000
## 4  0.8125000 0.8750000
## 5  0.9333333 0.9333333
## 6  0.8000000 0.8000000
## 7  0.6000000 0.5333333
## 8  0.9375000 0.9375000
## 9  0.7500000 0.9375000
## 10 0.8000000 0.7333333

     Pensando primeramente en poder aplicar un test \(\small{t}\) de Student para medidas pareadas, comenzamos comprobamos la normalidad de las diferencias entre las métricas:

# Test Shapiro-Wilk
shapiro.test(metricas$C50 - metricas$Rf)

## 
##  Shapiro-Wilk normality test
## 
## data:  metricas$C50 - metricas$Rf
## W = 0.83634, p-value = 0.03989

     En este caso, rechazamos con un 95% de confianza la hipótesis nula de que las diferencias siguen una distribución normal (\(\small{W}\) = 0,836, p-value = 0,039). Por tanto, debemos recurrir a una variante no paramétrica del test, por ejemplo el test de McNemar o el test de la suma de rangos con signo de Wilcoxon. Vamos a mostrar un ejemplo de aplicación de cada uno:

     Por una lado, el test de McNemar necesita de la matriz de confusión obtenida con las predicciones de los dos modelos considerados (comprobar que las parejas discordantes > 25), mientras que el test de Wilcoxon se aplica directamente sobre los datos de las métricas de validación cruzada:

# Matriz de confusion ejemplo
matriz_conf <- data.frame(si = c(29, 19), no = c(10, 137), 
                          row.names = c("si", "no"))

matriz_conf

##    si  no
## si 29  10
## no 19 137

# Test de McNemar
mcnemar.test(as.matrix(matriz_conf))

## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  as.matrix(matriz_conf)
## McNemar's chi-squared = 2.2069, df = 1, p-value = 0.1374

# Test de Wilcoxon
wilcox.test(metricas$C50, metricas$Rf, paired = TRUE)

## Warning in wilcox.test.default(metricas$C50, metricas$Rf, paired = TRUE): cannot
## compute exact p-value with ties

## Warning in wilcox.test.default(metricas$C50, metricas$Rf, paired = TRUE): cannot
## compute exact p-value with zeroes

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  metricas$C50 and metricas$Rf
## V = 7, p-value = 1
## alternative hypothesis: true location shift is not equal to 0

     Como podemos ver, con el test de Wilcoxon aparece un mensaje indicando que hay pliegues con empates y 0s. Nos basaremos, por tanto, en el resultado aportado por el test de McNemar, con el que no podemos evidenciar que existen diferencias significativas con un 95% de confianza entre el rendimiento de los dos modelos (\(\small{\chi^2_{Mc}}\) = 2,20 < \(\small{\chi^2_{1, 0.05}}\), p-value = 0,137).

DOS MODELOS EN DOMINIOS DISTINTOS

     En el caso de encontrarnos con la necesidad de comparar dos clasificadores sobre varios dominios o conjuntos de datos distintos. En cuanto a los test vistos anteriormente, el test \(\small{t}\) de Student asume que las medidas de rendimiento en los diferentes dominios a tratar tienen que ser comparables, y por otro lado el test de McNemar no es aplicable a más de dos clasificadores. La recomendación es aplicar el test de los rangos de signo de Wilcoxon.

Test de los rangos con signo de Wilcoxon

     Para muestras pareadas, el test de los rangos con signo de Wilcoxon es la alternativa no paramétrica del test \(\small{t}\) de Student cuando las muestras no siguen una distribución normal, por lo que compara las medianas en lugar de las medias.

HIPÓTESIS

\(\small{H_0}\): la mediana de las diferencias de cada par de medidas es igual a 0.

\(\small{H_1}\): la mediana de las diferencias de cada par de medidas es distinta de 0.

ESTADÍSTICOS

     Suponiendo que contamos con n dominios distintos, este test sirve para comparar diferentes pruebas realizadas (datasets) sobre la misma población con dos clasificadores distintos. El cálculo funciona de la siguiente manera:

• Siendo \(\small{p^i_A}\) y \(\small{p^i_B}\) las medidas de rendimiento de cada clasificador (clasificador A y clasificador B) en el dominio \(\small{i}\), se calculan las diferencias entre pares de medidas para cada dominio

• Se ordena el valor absoluto de di de menor a mayor y se le asigna un rango (±) a la diferencia.

• Se calculan dos sumas: la suma de los rangos con signo positivo (\(\small{T_+}\)) y la suma de los valores absolutos de los rangos con signo negativo (\(\small{T_-}\)). Las diferencias iguales a 0 se ignoran en estos cálculos). Si \(\small{H_0}\) fuera cierta, se esperaría que aproximadamente la mitad de las diferencias fueran positivas, y la otra mitad de diferencias negativas. Por tanto, la suma de rangos positivos sería aproximadamente igual a la suma de los rangos negativos (\(\small{T_+}\) ≈ \(\small{T_-}\)).

• Para examinar el supuesto anterior se calcula el estadístico (que sigue una distribución T de Wilcoxon) como el mínimo de las dos sumas de rangos

\({\small{T_{Wilcox} = min(T_+, T_-)}}\)

     Si n < 25, el p-value o la probabilidad de que el estadístico \({\small{T_{Wilcox}}}\) adquiera valores iguales o más extremos al observado se obtiene comparando el valor obtenido de \({\small{T_{Wilcox}}}\) con el valor crítico en una tabla Wilcoxon, para unos grados de libertad y significancia concretos. Si \({\small{T_{Wilcox}}}\) cae dentro del intervalo correspondiente, la diferencia no es significativa.

     En el caso de contar con más de 25 datos (n > 25), se puede asumir que el estadístico \({\small{T_{Wilcox}}}\) sigue una distribución aproximadamente normal, calculándose el siguiente estadístico:

donde \({\small{µ_{T_{Wilcox}}}}\) y \({\small{σ_{T_{Wilcox}}}}\) son la media y desviación estándar de la aproximación a la normal de la distribución \({\small{T_{Wilcox}}}\), en el caso de una H0 verdadera:

     En este caso, se rechaza \({\small{H_0}}\) si el \({\small{Z_{Wilcox}}}\) calculado es menor que el valor crítico de \({\small{Z}}\) de una tabla de distribución normal para unos grados de libertad y significancia concretos.

Adaptación a un único dominio

     El test de la suma de rangos con signo de Wilcoxon también podría aplicarse al caso de un único dominio (dataset) si se ha utilizado alguna técnica de remuestreo (validación cruzada, bootstraping, hold-out, etc.), pudiendo considerar el subconjunto generado en cada iteración como un dominio o dataset distinto.

Ejemplo en R

     Supongamos que contamos con las métricas de precisión de dos predictores sobre 10 conjuntos de datos de dominios distintos:

metricas <- data.frame(dataset = c("bladder", "breast", "colon", "kidney",
                                   "lung", "prostate", "thyroid", "oral",
                                   "uterine", "pancreatic"),
                       nnet = c(0.7541667, 0.8700000, 0.7150416, 0.8791667,
                                0.8361364, 0.6570000, 0.7133200, 0.8742437,
                                0.6834450, 0.6556061),
                       lda = c(0.7729167, 0.8566667, 0.7452914, 0.8697980,
                               0.8207008, 0.6660000, 0.7559288, 0.8714067,
                               0.6678606, 0.7420779))

metricas

##       dataset      nnet       lda
## 1     bladder 0.7541667 0.7729167
## 2      breast 0.8700000 0.8566667
## 3       colon 0.7150416 0.7452914
## 4      kidney 0.8791667 0.8697980
## 5        lung 0.8361364 0.8207008
## 6    prostate 0.6570000 0.6660000
## 7     thyroid 0.7133200 0.7559288
## 8        oral 0.8742437 0.8714067
## 9     uterine 0.6834450 0.6678606
## 10 pancreatic 0.6556061 0.7420779

Para compararlos, aplicaremos el test de la suma de rangos de Wilcoxon para muestras pareadas (ya que cada modelo se aplica sobre el mismo conjunto de datos):

# Test Wilcoxon
wilcox.test(metricas$nnet, metricas$lda, paired = TRUE)

## 
##  Wilcoxon signed rank exact test
## 
## data:  metricas$nnet and metricas$lda
## V = 19, p-value = 0.4316
## alternative hypothesis: true location shift is not equal to 0

Sobre los datos aplicados, el test no proporciona evidencias significativas con un 95% de confianza de una precisión distinta de los modelos aplicados (\({\small{W}}\) = 19 > \({\small{W^{10}_{0.05}}}\), p-value = 0,431).

VARIOS MODELOS EN DOMINIOS DISTINTOS

     A la hora de comparar múltiples modelos, podríamos caer en la tentación de repetir, por ejemplo, un test \({\small{t}}\) de Student para cada par, sin embargo, esto supondría un aumento de la probabilidad de cometer un Error de Tipo I, sobre todo para muchas comparaciones. En su lugar, se aplican test estadísticos (paramétricos o no) que realizan contrastes de varias hipótesis al mismo tiempo. Con ellos podemos detectar si al menos el rendimiento de dos modelos son significativamente diferentes (pero no nos aportan información directa sobre cuales).

ANOVA de una vía para medidas repetidas

     Como extensión del test \({\small{t}}\) de Student para medidas pareadas, el test ANOVA de medidas repetidas (one way repeated measures ANOVA) nos permite comparar las diferencias entre las medias de las medidas de rendimiento de más de dos modelos entre los diferentes dominios o grupos de datos, y determinar si estas diferencias son estadísticamente significativas. Estas medias las compara mediante el estudio de las varianzas, concretamente, se calcula la media de cada grupo (modelo) para comparar la varianza de estas medias (intervarianza) frente a la varianza promedio dentro de cada modelo (intravarianza).

HIPÓTESIS

\(\small{H_0}\): las medias del rendimiento de todos los modelos son iguales (\(\small{µ_0}\) = \(\small{µ_1}\) = … = \(\small{µ_n}\))

\(\small{H_1}\): al menos dos medias son significativamente distintas

CONDICIONES

• Normalidad: las muestras deben ser independientes y proceder de una distribución aproximadamente normal.

• Esfericidad: las varianzas de las diferencias entre todos los pares de medidas repetidas deben ser similares (equivalente a la homocedasticidad para el ANOVA de medidas pareadas). La esfericidad puede comprobarse aplicando el test de Mauchly.

     Tener en cuenta que los conjuntos de datos deberían tener aproximadamente el mismo tamaño y que las medidas de rendimiento de cada modelo deben tener la misma escala.

ESTADÍSTICO

     Un ANOVA de medidas repetidas calcula un estadístico \(\small{F}\) de una forma similar a un ANOVA de medidas independientes. Entre los cálculos implicados para la obtención del estadístico, encontramos sumas de cuadrados, grados de libertad (df) y cuadrados medios (MS):

• \(\small{SS_{total}}\): suma de cuadrados total (total de la varianza en los datos). Se corresponde con la suma de los cuadrados de las diferencias de cada medida respecto a la media global de todos los datos. Los grados de libertad corresponden al número de medidas (\(\small{n}\)) por número de grupos o modelos (\(\small{k}\)), menos 1: (\(\small{n}\) x \(\small{k}\)) -1.

• \(\small{SS_{intra}}\): suma de cuadrados dentro de cada grupo. Grados de libertad = \(\small{n}\) · (\(\small{k - 1}\)).

• \(\small{SS_{modelo}}\): varianza explicada por el modelo. Grados de libertad = \(\small{k - 1}\).

• \(\small{SS_{error}}\): varianza no explicada por el modelo (\(\small{SS_{intra}}\) – \(\small{SS_{modelo}}\)). Grados de libertad = (\(\small{n - 1}\)) · (\(\small{k - 1}\)).

• \(\small{df_{modelo}}\) = \(\small{k - 1}\)

• \(\small{df_{error}}\) = (\(\small{k - 1}\)) · (\(\small{n - 1}\))

• \(\small{MS_{modelo}}\) = \(\small{SS_{modelo}}\) / \(\small{df_{modelo}}\)

• \(\small{MS_{error}}\) = \(\small{SS_{error}}\) / \(\small{df_{error}}\)

• \(\small{F}\) = \(\small{MS_{modelo}}\) / \(\small{MS_{error}}\)

     El estadístico \(\small{F}\) sigue una distribución F Fisher-Snedecor. Si la probabilidad de obtener valores del estadístico como el obtenido o mayores es menor al α establecido, se rechaza \(\small{H_{0}}\).

     Al igual que otros test vistos, el resultado del ANOVA nos indica si existen diferencias entre algún par de modelos, pero no nos indica explícitamente cuales.

Test de Friedman

     Aunque el test ANOVA para medidas repetidas resulta robusto frente a la violación de la condición de normalidad, en ocasiones la esfericidad puede ser difícil de comprobar para comparar clasificadores, o puede que contemos con medidas de rendimiento categóricas. Por tanto, como alternativa no paramétrica, podemos optar por el test de Friedman, el cual se considera una extensión de la prueba de los rangos con signo de Wilcoxon para más de dos grupos. La comparación será entre las medianas en lugar de las medias.

     Incluso cumpliéndose las condiciones, puede darse el caso de que los resultados del ANOVA y el test de Friedman den resultados similares.

HIPÓTESIS

\(\small{H_0}\): todas las medianas son iguales

\(\small{H_1}\): al menos dos medianas difieren

Nota: Al igual que ocurre con el test de Wilcoxon, el test de Friedman se basa en los rangos de cada clasificador más que en las medidas de rendimiento.

ESTADÍSTICO

     Se calcula también un estadístico \(\small{F_r}\) con distribución \(\small{\chi^2}\). El procedimiento a groso modo es el siguiente:

1. Dado un conjunto de \(\small{k}\) grupos (modelos) y \(\small{n}\) bloques (medidas), se obtiene el rango de los datos (medida de rendimiento) por cada bloque. Es decir, se asignan los rangos por filas de datos.
2. Por cada grupo o columna se suman los rangos (\(\small{T_1, T_2, ..., T_k}\)).
3. Se calcula el estadístico \(\small{F_r}\) como

4. Obtener el p-value para \(\small{k - 1}\) grados de libertad.

Comparaciones post-hoc

     Los test comentados solo nos aportan información sobre la existencia de diferencias significativas, pero no entre qué modelos se dan estas diferencias. Para determinar esto último, hemos de recurrir a test post-hoc, de los cuales también existen versiones paramétricas y no paramétricas.

     En caso de que el ANOVA de medidas repetidas sea significativo, aplicaríamos uno de los siguientes:

PARAMÉTRICOS

• Test de Tukey: tiene menos probabilidad de cometer un Error de Tipo I que el test \(\small{t}\).

• Test de Dunnet: cuando las comparaciones no son dos a dos, sino de todos los clasificadores frente a uno control.

• Test de Bonferroni: para un número de comparaciones pequeño, de lo contrario tiende a ser conservador.

• Test de Bonferroni-Dunn: divide α por el número de comparaciones a realizar. Corrige el conservadurismo del test de Bonferroni.

     En el caso de que el test de Friedman proporcione resultados significativos, aplicaríamos uno de los siguientes:

NO PARAMÉTRICOS

• Test de Nemenyi

• Otros: test de Holm, test de Hommel, test de Hochberg

Ejemplo en R

Caso paramétrico

     Para este ejemplo supondremos que contamos con las métricas de precisión de 3 modelos sobre los 10 conjuntos de datos de distintas fuentes:

metricas <- data.frame(nnet = c(0.95343, 0.95428, 0.95637, 0.95334, 0.95795,
                                0.95935, 0.95421, 0.95617, 0.95462, 0.95640),
                       lda = c(0.85374, 0.82696, 0.78597, 0.83563, 0.84223,
                               0.75414, 0.85768, 0.84611, 0.80011, 0.83464),
                       svm = c(0.93707, 0.95418, 0.93346, 0.94629, 0.95895,
                               0.90717, 0.95892, 0.96263, 0.97544, 0.96700),
                       row.names = c("bladder", "breast", "colon", "kidney",
                                   "lung", "prostate", "thyroid", "oral",
                                   "uterine", "pancreatic"))

metricas

##               nnet     lda     svm
## bladder    0.95343 0.85374 0.93707
## breast     0.95428 0.82696 0.95418
## colon      0.95637 0.78597 0.93346
## kidney     0.95334 0.83563 0.94629
## lung       0.95795 0.84223 0.95895
## prostate   0.95935 0.75414 0.90717
## thyroid    0.95421 0.85768 0.95892
## oral       0.95617 0.84611 0.96263
## uterine    0.95462 0.80011 0.97544
## pancreatic 0.95640 0.83464 0.96700

# Cambio de formato wide a long: concatemanos los valores de las metricas de cada modelo en un unico vector, con un segundo vector indicando de qué modelo
# procede el dato y un tercero que indique el dataset
metricas_stack <- stack(metricas)
metricas_stack$dataset <- as.factor(rep(row.names(metricas), times = 3))
names(metricas_stack) <- c("Accuracy","Model", "Dataset")
head(metricas_stack)

##   Accuracy Model  Dataset
## 1  0.95343  nnet  bladder
## 2  0.95428  nnet   breast
## 3  0.95637  nnet    colon
## 4  0.95334  nnet   kidney
## 5  0.95795  nnet     lung
## 6  0.95935  nnet prostate

     La opción paramétrica con la que contamos para este caso es el test ANOVA de una vía para medidas pareadas. Comenzamos comprobando las condiciones para su uso (normalidad de las métricas de cada modelo en cada conjunto de datos, y la esfericidad de la matriz de covarianzas):

# Test Saphiro Wilk (comprobar normalidad)
tapply(metricas_stack$Accuracy, metricas_stack$Model, shapiro.test)

## $nnet
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92026, p-value = 0.3591
## 
## 
## $lda
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.87423, p-value = 0.1119
## 
## 
## $svm
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92585, p-value = 0.4084

     La normalidad no puede ser rechazada en ninguno de los casos según el test de Saphiro Wilk con un 95% de confianza (\(\small{W_{nnet}}\) = 0,92, p-value = 0,359; \(\small{W_{lda}}\) = 0,87, p-value = 0,111; \(\small{W_{svm}}\) = 0,92, p-value = 0,408).

library(car)

# Test de Mauchly (comprobar esfericidad). Se utiliza la funcion Anova del 
# paquete "car", que devuelve dicho test:

# factors defining the intra-subject model for multivariate repeated-measures
# data
modelos <- factor(c("nnet", "lda", "svm"))

options(contrasts = c("contr.sum", "contr.poly"))

# Anova
summary(Anova(lm(as.matrix(metricas) ~ 1), # modelo lineal
              idata = data.frame(modelos), 
              idesign = ~ modelos),
        multivariate = FALSE)

## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##              Sum Sq num Df  Error SS den Df  F value    Pr(>F)    
## (Intercept) 24.8329      1 0.0065562      9 34089.25 < 2.2e-16 ***
## modelos      0.1113      2 0.0070349     18   142.43 9.252e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##         Test statistic p-value
## modelos        0.67869 0.21217
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##          GG eps Pr(>F[GG])    
## modelos 0.75682  2.299e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            HF eps Pr(>F[HF])
## modelos 0.8773595 1.4889e-10

     No existen evidencias significativas para rechazar la esfericidad de la matriz de covarianzas (p-value = 0,212).

     Ante el cumplimiento de la normalidad y la esfericidad, podemos aceptar los resultados del ANOVA para medidas repetidas, el cual indica la existencia de diferencias significativas entre los tres modelos (p-value = \(\small{9,25x10^{-12}}\)), pero no nos indica cuales difieren. Para ello optaremos por aplicar el test de post hoc de Tukey. Sin embargo, con objetos generados por una función ANOVA para medidas repetidas no podemos aplicar las funciones de R disponibles para la aplicación directa del test de Tukey. Una forma de hacerlo es utilizando la función glht() para hipótesis lineales generales y comparaciones múltiples sobre modelos paramétricos del paquete multcomp:

library(nlme)
# Modelo lineal generalizado que compensa el hecho de que no puede existir 
# variabilidad entre los grupos
modelo_lme <- lme(Accuracy ~ Model, random = ~1 | Dataset, 
                  data = metricas_stack)

library(multcomp)
# Test Tukey, aplicado al modelo lineal generalizado
tukey_lme <- glht(modelo_lme, linfct=mcp(Model = "Tukey"))
summary(tukey_lme)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lme.formula(fixed = Accuracy ~ Model, data = metricas_stack, 
##     random = ~1 | Dataset)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## lda - nnet == 0 -0.131891   0.008841 -14.918   <1e-04 ***
## svm - nnet == 0 -0.005501   0.008841  -0.622    0.808    
## svm - lda == 0   0.126390   0.008841  14.296   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

El test muestra diferencias altamente significativas entre los modelos LDA-NNET (\(\small{Z}\) = -14,91, p-value = \(\small{1x10^{-5}}\), siendo NNET mejor que LDA ya que el estimador de la media de las diferencias es negativo), y entre SVM-LDA (\(\small{Z}\) = 14,29, p-value = \(\small{1x10^{-5}}\), siendo SVM mejor que LDA ya que el estimador de la media de las diferencias es positivo). Con esto podemos concluir que el modelo LDA es el peor de los tres.

Caso no paramétrico

     Supongamos que contamos con un caso igual al anterior, habiendo obtenido unas métricas distintas:

metricas <- data.frame(nnet = c(0.3669627, 0.8487879, 0.2656273, 0.8660868,
                                0.7509778, 0.2667577, 0.4762539, 0.8290359,
                                0.3857618, 0.3987424),
                       lda = c(0.4325550, 0.8287879, 0.3012339, 0.8553988,
                               0.7186083, 0.2873961, 0.5731759, 0.8234315,
                               0.3661388, 0.5699495),
                       svm = c(0.3498855, 0.8587879, 0.3365727, 0.8917666,
                               0.7527696, 0.2972005, 0.4869577, 0.8400313,
                               0.3945373, 0.6073158),
                       row.names = c("bladder", "breast", "colon", "kidney",
                                   "lung", "prostate", "thyroid", "oral",
                                   "uterine", "pancreatic"))
metricas

##                 nnet       lda       svm
## bladder    0.3669627 0.4325550 0.3498855
## breast     0.8487879 0.8287879 0.8587879
## colon      0.2656273 0.3012339 0.3365727
## kidney     0.8660868 0.8553988 0.8917666
## lung       0.7509778 0.7186083 0.7527696
## prostate   0.2667577 0.2873961 0.2972005
## thyroid    0.4762539 0.5731759 0.4869577
## oral       0.8290359 0.8234315 0.8400313
## uterine    0.3857618 0.3661388 0.3945373
## pancreatic 0.3987424 0.5699495 0.6073158

# Cambio de formato wide a long: concatemanos los valores de las metricas de cada modelo en un unico vector, con un segundo vector indicando de qué modelo
# procede el dato y un tercero que indique el dataset
metricas_stack <- stack(metricas)
metricas_stack$dataset <- as.factor(rep(row.names(metricas), times = 3))
names(metricas_stack) <- c("Accuracy","Model", "Dataset")
head(metricas_stack)

##    Accuracy Model  Dataset
## 1 0.3669627  nnet  bladder
## 2 0.8487879  nnet   breast
## 3 0.2656273  nnet    colon
## 4 0.8660868  nnet   kidney
## 5 0.7509778  nnet     lung
## 6 0.2667577  nnet prostate

     Para saber qué tipo de test podemos aplicar, comenzamos comprobando la normalidad de de las métricas de cada modelo en cada conjunto de datos:

# Test Saphiro Wilk
tapply(metricas_stack$Accuracy, metricas_stack$Model, shapiro.test)

## $nnet
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.83722, p-value = 0.04086
## 
## 
## $lda
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.89558, p-value = 0.1958
## 
## 
## $svm
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.87459, p-value = 0.113

     Según el test, las métricas obtenidas por el modelo NNET no siguen una distribución normal (\(\small{W_{nnet}}\) = 0,837, p-value = 0,04). Sin necesidad de comprobar la esfericidad, procederemos a aplicar la alternativa no paramétrica al ANOVA de una vía para medidas repetidas, el test de Friedman:

# Test de Friedman
friedman.test(Accuracy ~ Model | Dataset, data = metricas_stack)

## 
##  Friedman rank sum test
## 
## data:  Accuracy and Model and Dataset
## Friedman chi-squared = 7.4, df = 2, p-value = 0.02472

     Esta vez el test ha encontrado evidencias significativas (\(\small{F}\) = 7,4, p-value = 0,024), con un 95% de confianza, de que existen diferencias entre la precisión de los tres modelos con los 10 conjuntos de datos utilizados. Para saber entre qué modelos existen estas diferencias, aplicaremos el test post hoc de Nemenyi:

library(PMCMR)

# Test Nemenyi
posthoc.friedman.nemenyi.test(Accuracy ~ Model | Dataset, data = metricas_stack)

## 
##  Pairwise comparisons using Nemenyi multiple comparison test 
##              with q approximation for unreplicated blocked data 
## 
## data:  Accuracy and Model and Dataset 
## 
##     nnet  lda  
## lda 0.973 -    
## svm 0.037 0.065
## 
## P value adjustment method: none

     El test nos devuelve los p-values de las comparaciones entre pares de modelos. En este caso, no se encuentran diferencias significativas entre SVM-LDA (p-value = 0,065), ni entre LDA-NNET (p-value = 0,973), pero sí entre SVM-NNET (p-value = 0,037), teniendo SVM una precisión ligeramente superior.

BIBLIOGRAFÍA

Brian S. Everitt and Torsten Hothorn. A Handbook of Statistical Analyses Using R, 2009.

Raschka S. Model Evaluation, Model Selection, and Algorithm Selection in Machine Learning. University of Wisconsin-Madison, Department of Statistics (2018).

Taskin Kavzoglu. Handbook of Neural Computation, 2017. Robert H. Riffenburgh, Statistics in Medicine (Second Edition), 2006.

This work by Cristina Gil Martínez is licensed under a Creative Commons Attribution 4.0 International License.