Este guion está diseñado para que vaya aplicando los test de hipótesis que se han visto en la parte teórica del curso. Su realización se hará durante las clases presenciales del Máster. El alumno no tiene que entregar nada relativo a esta parte. De hecho, el 95 por ciento de este guion está resuelto.

Se usarán varios conjuntos de datos cuyos enlaces se proporcionan a lo largo del guion conforme se van necesitando. También puede descargarlos todos desde el siguiente enlace:

https://decsai.ugr.es/jccubero/Master/IntrodEstadisticaData.zip

En primer lugar se verá el uso de la herramienta Jamovi que es un entorno gráfico para seleccionar los tests de hipótesis que se quieran ejecutar. Si el alumno termina durante las clases prácticas esta parte, puede pasar a la parte de R en la que se ve cómo ejecutar los mismos tests (y otros más avanzados) desde R.

En lo que sigue, nos vamos a centrar en el caso de una variable dependiente continua (puede ser ordinal en el caso de los tests no paramétricos) y una única variable independiente o factor. En el caso de que tuviésemos en cuenta varios factores, habría que tener en cuenta la interacción entre los distintos niveles o grupos de dichas variables (multiway Anova). Si tuviésemos varias variables dependientes tendríamos que recurrir a las técnicas englobadas bajo el término MANOVA. Si consideramos también variables continuas como factores, recurriríamos al ANCOVA.

Resumen de los tests que se van a ver en este guion (con Jamovi y con R):

Tests para muestras independientes

Las filas han de ser independientes entre sí. Los grupos han de ser independientes entre sí.
- Test para dos muestras
  - Paramétrico. Los residuos han de ser normales. Dependiendo de si las varianzas son similares o no (según nos indique el test de Levene) tenemos:
    - Varianzas similares: t-test
    - Varianzas distintas: Welch
  - No paramétrico: Mann-Whitney (es un test sobre las medianas cuando las distribuciones son similares y es un test de dominancia estocástica en otro caso)
- Test para k muestras
  - Paramétrico. Los residuos han de ser normales.
    - Varianzas similares: One-way Anova. Post-hoc: Tukey HSD
    - Varianzas distintas: Welch. Post-hoc: Games-Howell
    - En cualquiera de los casos anteriores se puede aplicar un post-hoc genérico como por ejemplo Holm.
  - No paramétrico: Kruskal-Wallis (es un test sobre las medianas cuando las distribuciones son similares y es un test de dominancia estocástica en otro caso). Post-hoc: Dwass-Steel-Critchlow-Fligner, Dunn, Conover, Nemenyi. También se puede aplicar un post-hoc genérico.
Tests para muestras dependientes (medidas repetidas)

Las filas han de ser independientes entre sí.
- Test para dos muestras
  - Paramétrico. Los residuos de las diferencias han de ser normales. Se aplica el paired t-test
  - No paramétrico: Si la diferencia de las dos variables es una distribución simétrica, se aplicará Wilcoxon. En otro caso, el sign-test
- Test para k muestras
  - Paramétrico. Los residuos han de ser normales. Se aplica el Anova de medidas repetidas. Si no se verifica la hipótesis de esfericidad (test de Mauchly) se aplica la corrección de GreenHouse-Geiser. Post-hoc con Tukey o uno genérico.
  - No paramétrico: Friedman (es un test sobre las medianas cuando las distribuciones son similares y es un test de dominancia estocástica en otro caso). Post-hoc: Durbin-Conover, Nemenyi, Eisinga, etc. También se puede aplicar un post-hoc genérico.

Puede consultar en el siguiente enlace una tabla completa con los distintos tests a aplicar, dependiendo de los datos con los que trabajamos:

http://statkat.com/stattest_overview.php

1 Tests de Hipótesis con Jamovi

Puede descargar Jamovi desde el siguiente enlace:

https://www.jamovi.org/

1.1 Tests para muestras independientes

1.1.1 Prueba para dos muestras independientes (test paramétrico)

dataset: Datos de empleados
test: ¿Puede considerarse que el salario medio de los hombres es igual al de las mujeres?
Variable dependiente: salario
Variable independiente (o factor): sexo

El conjunto de datos de empleados es un dataset ficticio incluido en la distribución de SPSS y almacena información de los trabajadores de una empresa.

Tenemos un conjunto de datos de hombres (muestra 1) y un conjunto de datos de mujeres (muestra 2). Para obtener dichas muestras basta seleccionar como variable dependiente el salario y como variable de agrupación el sexo. Lo hacemos desde “ANOVA / One-Way Anova”. Vamos a aplicar un test de hipótesis para comparar ambas muestras. Observe que para aplicar el test no es necesario que ambas muestras sean del mismo tamaño (de hecho, no hay el mismo número de hombres que de mujeres) En una primera instancia, aplicamos el test basado en la t-Student (T-test), aunque posteriormente veremos que algunas condiciones que exige este test no se verifican y tendremos que aplicar otros tests alternativos.

El test es significativo por lo que rechazaríamos la hipótesis nula y aceptaríamos que los hombres no tienen el mismo salario que las mujeres. En el gráfico descriptivo, podemos apreciar la diferencia que hay entre las medias de los dos grupos (la barra sobre el valor de la media representa el intervalo de confianza para dicho valor)

En cualquier caso, tenemos que analizar los requisitos del test:

Independencia de las dos muestras: En el diseño de un experimento, el científico debe asegurar que los sujetos (personas, animales, materiales, etc.) son elegidos y asignados a cada muestra de forma aleatoria. Esto es lo que garantiza que los sujetos de cada muestra sean independientes y los resultados del experimento no se vean condicionados por un posible sesgo en la elección de los sujetos o la asignación de un sujeto a una muestra u otra. Por ejemplo, podríamos estar interesados en ver el efecto de un relajante en una población de ratones, para lo cual medimos el número de horas de sueño. Normalmente se usarán dos muestras: una de ellas recibirá el relajante y la otra no (será un grupo de control). En primer lugar, los ratones deberían elegirse de forma aleatoria, de forma que no deberían seleccionarse únicamente los machos, o los que tengan una característica determinada. En segundo lugar, también ha de ser aleatoria la asignación de cada ratón al grupo de los que recibirán el relajante o al grupo de control. En definitiva, se ha de asegurar:
- La independencia entre las filas de un mismo grupo
- La independencia entre las filas de distintos grupos
En nuestro caso, no vamos a extraer una muestra de una población de posibles sujetos sino que vamos a trabajar con todos los datos que ya han sido obtenidos previamente. En cualquier caso, debemos asegurar que los sujetos que vamos a comparar (las filas o registros del dataset) son independientes. Supongamos por ejemplo los datos de los empleados.
- ¿Podemos garantizar que la elección de cada sujeto ha sido aleatoria? Debemos suponer que los empleados actuales de la empresa representan una muestra de un conjunto mucho mayor de posibles empleados que, cumpliendo los requisitos, podrían haber optado a un puesto. ¿La incorporación a la empresa de un empleado viene condicionada por la contratación de otro? Si fuese así, no podríamos considerar que las filas son independientes entre sí.
- ¿Podemos garantizar que la muestra constituida por los empleados varones es independiente de la muestra de las mujeres? Supongamos que la mayor parte de los empleados de la empresa pertenecen a grupos familiares y que al contratar a una persona, se intenta siempre contratar también a un hermano de otro sexo de la misma familia. En este caso, es muy posible que exista un sesgo determinado por el grupo familiar. Por ejemplo, puede ser que los miembros de una misma familia tengan un nivel de estudios similar y por tanto es muy posible que tengan la misma categoría laboral, salarios similares, etc. Esta dependencia introduce un sesgo que puede invalidar el resultado del test ya que podría ser no significativo (no se rechaza que las medias son similares) cuando realmente, si eliminamos los datos de los hermanos sí obtendríamos un test significativo (se rechararía por tanto que las medias son similares).
En todos los tests de hipótesis que vamos a ver en este apartado de muestras independientes supondremos que se verifican ambas hipótesis de independencia y por tanto, la presencia de una fila no condiciona la presencia de otra fila ni en el mismo grupo ni entre grupos diferentes.

Es importante destacar que hemos denominado a salario como variable dependiente. Esto no contradice lo dicho anteriormente. Nuestro objetivo es ver si las medias del salario son las mismas en las dos muestras (hombre y mujeres) Si rechazamos dicha hipótesis, estaremos admitiendo que el sexo influye a la hora de determinar el salario: el salario depende del sexo. Para llegar a esta conclusión necesitamos que la muestra de los hombres y de las mujeres no haya sido obtenida aplicando algún tipo de sesgo, es decir, deben ser independientes.
Homocedasticidad (homogeneidad de las varianzas): Todos los grupos han de tener una variabilidad similar en cuanto a la variable dependiente. El gráfico de los diagramas de cajas nos dice que no se cumple el requisito de homogeneidad de varianzas ya que la altura de la caja de los hombres es bastante mayor que la de las mujeres. El test de Levene (homogeneity) confirma lo anterior ya que rechaza la hipótesis nula de igualdad de varianzas.

En el caso de varianzas distintas, aplicamos mejor el test de Welch, seleccionando la casilla correspondiente:

El test resulta significativo. En cualquier caso, el test de Welch también requiere que se cumpla la hipótesis de normalidad:

Normalidad de los residuos. Si seleccionamos las casillas de “Normality test” y “q-q plot” en el apartado “Assumption checks” obtenemos el gráfico que aparece abajo. Podemos apreciar que los residuos NO siguen una distribución normal (los puntos deberían estar todos cercanos a la línea de referencia) El correspondiente test de Shapiro Wilks también rechaza la hipótesis de normalidad.

En resumen: Hemos aplicado el test de Welch al tener los grupos varianzas distintas. Éste resulta significativo por lo que podemos concluir que, efectivamente, el salario medio de los hombres es distinto al de las mujeres. Como la media del salario de los hombres es 41442 y el de las mujeres es 26032, podemos concluir que el salario medio de los hombres es mayor estricto que el de las mujeres. Esta afirmación no tiene una garantía estadística completa ya que, aunque los tests paramétricos como la t-Student y el de Welch son bastante robustos a la violación del requisito de la normalidad, es conveniente aplicar un test no paramétrico que, aunque tenga menos potencia, no requiere que se cumpla la hipótesis de normalidad. Esto lo haremos a continuación.

Nota: Estamos planteando como hipótesis nula la igualdad de las medias y como alternativa que sean distintas (test a dos colas). Es posible plantear tests más específicos en los que la hipótesis nula sea, por ejemplo, que la media de los hombres es mayor y la alternativa sea que la media sea menor o igual (test a una cola)

1.1.2 Test no paramétrico para comparar dos muestras independientes

dataset: Datos de empleados
test (en el caso de varianzas similares): ¿Puede considerarse que la mediana del salario de los hombres es igual al de las mujeres?
test (en el caso de varianzas distintas): ¿Puede considerarse que el salario de los hombres es estocásticamente igual al de las mujeres?
Variable dependiente: salario
Variable independiente (o factor): sexo

En el mismo cuadro de diálogo del apartado paramétrico , seleccionamos la casilla Mann-Whitney y nos debe salir el siguiente resultado:

El test de Mann-Whitney resulta significativo. Veamos los requisitos del test:

Como en todos los tests de este guion, asumimos que las filas son independientes entre sí.
El test de Mann-Whitney no impone el requisito de que las distribuciones de los grupos sean normales. En cualquier caso, si queremos que sea un test sobre las medianas, sí se necesita que las distribuciones de la variable dependiente en cada uno de los grupos sean similares (el histograma del salario en los hombres debe ser similar al histograma del salario en las mujeres) Pero resulta que ni siquiera las varianzas de los grupos son similares, por lo que no podemos asegurar que se verifica este requisito. Por lo tanto, el test de Mann-Whitney no podemos decir que sea un test sobre las medianas sino un test de dominancia estocástica.

En resumen, podemos concluir que el salario de los hombres es estocásticamente mayor que el de las mujeres, es decir, que si elegimos al azar un hombre y una mujer, es más probable que el salario mayor corresponda al del hombre. No podemos concluir que las medianas sean distintas.

1.1.3 Anova para k muestras independientes (test paramétrico)

1.1.3.1 Test conjunto (omnibus)

dataset: Datos de empleados
test: ¿Puede considerarse que el salario medio es el mismo en cada una de las tres categorías laborales (k=3)?
Variable dependiente: salario
Variable independiente (o factor): catlab

Seleccionamos ANOVA/One-Way Anova

El test resulta significativo, pero debemos comprobar si se verifican los requisitos del test que, al igual que ocurría en el caso de dos muestras independientes, son la independencia de los datos, la homocedasticidad y la normalidad.

Independencia de los datos. Tal y como dijimos al inicio de la sección, suponemos que todas las filas son independientes dentro de cada grupo y entre grupos distintos.
Homocedasticidad. Podemos apreciar en el gráfico de cajas que es muy distinta en los tres grupos de la categoría laboral. Claro está, el test de Levene rechaza la igualdad de varianzas. Por eso hemos aplicado el test de Welch en vez de un Anova asumiendo varianzas similares (Fisher)

Normalidad de los residuos. Si analizamos los residuos, vemos que una parte de éstos se alejan mucho de la línea de referencia por lo que no podemos asumir que éstos sigan una distribución normal.

En resumen, al tener varianzas distintas en cada uno de los grupos, hemos aplicado el test Anova de Welch y ha resultado significativo, por lo que podemos concluir que hay al menos un grupo con un salario medio distinto al resto de los grupos. En cualquier caso, aunque tal y como ocurría en el caso de dos muestras independientes el test de Anova es bastante robusto a violaciones de la normalidad de los residuos, este requisito no se cumple por lo que no tenemos una garantía estadística completa y por tanto es preferible ejecutar el mismo estudio pero usando un test no paramétrico (lo haremos posteriormente)

1.1.3.2 Post hoc

Procedemos a realizar las comparaciones múltiples dos a dos entre todos los grupos de la variable independiente o factor (categoría laboral) Si tenemos \(k\) grupos, en el post hoc estaremos planteando un total de \(k(k-1)/2\) comparaciones, por lo que, tal y como hemos visto en el curso, los métodos tendrán que controlar el error FWER. Los métodos que hemos denominado específicos realizan algún tipo de penalización interna en la construcción del propio test y otros (los que hemos denominado genéricos) aplican un test normal para el caso de dos muestras a cada una de las \(k(k-1)/2\) posibles combinaciones y luego aplican una penalización del nivel de significación (Bonferroni, Sidak, Holm, Hochberg, etc.)

Uno de los métodos específicos más usados en el caso del Anova paramétrico con varianzas distintas es el test de Games-Howell. Si las varianzas no fuesen distintas, seleccionaríamos el método de Tukey HSD (aunque el test de Games-Howell también tiene una potencia similar a éste en el caso de varianzas iguales)

Nos debe salir lo siguiente:

Vemos que las comparaciones de todos los pares de categoría (Administrativo vs Seguridad, Administrativo vs Directivo, Seguridad vs Directivo) han resultado significativas por lo que concluimos que las tres categorías laborales tienen medias distintas.

Para que compruebe la importancia de considerar la varianza, ejecute el mismo post hoc pero usando el método de Tukey. Puede observar que este test no es capaz de rechazar la igualdad entre las medias de Administrativo y Seguridad.

Tal y como hemos comentado anteriormente, para realizar cada una de las \(k(k-1)/2\) comparaciones múltiples, también podríamos haber aplicado un test genérico sobre cada par de medias y penalizar los niveles de significación con algún método como Holm. Para lanzar este tipo de tests, nos tenemos que ir a otro cuadro de diálogo en Jamovi. Debemos seleccionar ANOVA/ANOVA y en la pestaña Post Hoc Tests seleccionamos catlab. Usamos el método de penalización por pasos de Holm:

Nos debe salir lo siguiente:

Vemos que el test sí distingue entre directivos y las otras dos categorías pero no es capaz de rechazar la igualdad entre administrativos y seguridad, algo que sí hacía el test de Games-Howell. La razón es que Jamovi lanza el post hoc usando como test genérico para realizar cada comparación el t-test. Pero para estimar la varianza utiliza la denominada varianza agregada (pooled variance). Esto sólo tiene sentido cuando las varianzas de cada grupo son similares, algo que ya hemos visto que no ocurre. Por lo tanto, en los cómputos del estadístico, los directivos aportan una varianza muy alta que domina sobre el resto de los grupos. En definitiva, en este caso de varianzas distintas no deberíamos aplicar este post hoc y deberíamos aplicar únicamente el test de Games-Howell. Cuando veamos la parte del guion con R, veremos cómo aplicar tests genéricos de forma correcta.

En resumen, si va a usar Jamovi, cuando los datos no cumplan el requisito de homocedasticidad, use únicamente el test de Games-Howell como post-hoc al Anova de Welch

Nota El test de Welch sólo está disponible en el cuadro de díalogo ANOVA/One-way Anova pero no lo está bajo el cuadro de diálogo de ANOVA/ANOVA. La razón es que el ANOVA genérico permite utilizar varias variables independientes o factores para analizar cómo interactúan entre ellas (el Anova de una vía -one way Anova- sólo considera un único factor) En estos casos, el test de Welch no es aplicable.

1.1.4 Anova para k muestras independientes (test no paramétrico)

1.1.4.1 Test conjunto (omnibus)

dataset: Datos de empleados
test (en el caso de varianzas distintas): ¿Puede considerarse que el salario en cada una de las tres categorías laborales es estocásticamente igual al del resto de categorías?
test (en el caso de varianzas similares): ¿Puede considerarse que la mediana del salario es la misma en cada una de las tres categorías laborales?

En el Anova paramétrico vimos que no se cumplía el requisito de normalidad de los residuos. Si bien es cierto que el Anova es bastante robusto a desviaciones de la normalidad, si queremos guardarnos las espaldas, podemos lanzar un test no paramétrico que no imponga dicho requisito. El test no paramétrico a ejecutar sería el de Kruskal-Wallis.

Nos debe salir lo siguiente:

Al igual que en el test no paramétrico de 2 muestras independientes, si las varianzas (en realidad el histograma completo) de cada grupo fuesen similares, el test no paramétrico de Kruskal-Wallis sería un test sobre las medianas por lo que podríamos concluir que hay alguna mediana significativamente mayor (o menor). Como las varianzas son realmente muy distintas en cada grupo, el test sólo puede concluir que hay algún grupo que domina estocásticamente sobre otro.

1.1.4.2 Post hoc

Si en el mismo cuadro de diálogo del test de Kruskal-Wallis seleccionamos el post-hoc de DSCF (Dwass-Steel-Critchlow-Fligner), obtenemos:

por lo que podemos concluir que el grupo de los directivos domina estocásticamente sobre los otros dos y los mismo ocurre con el grupo de administrativos frente a los de seguridad. Así pues, por ejemplo, para este último caso, si elegimos al azar un empleado administrativo y un empleado de seguridad, es más probable que el salario menor corresponda al de administrativo. Recordemos que si los grupos hubiesen tenido distribuciones similares, también podríamos concluir que las medianas son distintas.

1.2 Tests para muestras dependientes o de medidas repetidas

Ahora vamos a trabajar con muestras dependientes. El ejemplo típico es cuando a un mismo sujeto se le realizan varias medidas (de ahí el nombre de medidas repetidas) Por ejemplo, se le administra un fármaco y se mide alguna variable que nos diga cómo responde al tratamiento. Si realizamos tres mediciones (separadas entre sí un número de días), tendremos tres columnas de datos. Las columnas representan muestras dependientes ya que el valor de un registro en cualquiera de ellas está obviamente condicionado por el valor del mismo registro en las otras columnas. En cualquier caso, debemos seguir exigiendo que la elección de los sujetos haya sido aleatoria, es decir, las filas han de ser independientes.

En resumen, en esta sección vamos a aplicar tests en los que:

Las filas dentro de un grupo son independientes entre sí
Las filas de un grupo no son independientes del resto de los grupos.

Al igual que en el caso de muestras independientes, tendremos tests paramétricos que exigen ciertos requisitos y tests no paramétricos que rebajan la exigencia de algunos.

1.2.1 Prueba T (t-Student) para dos muestras pareadas (medidas repetidas)

dataset: Datos de empleados
test: ¿Puede considerarse que la diferencia entre las medias del salario actual y el salario inicial es cero?

Este es un ejemplo claro de variables pareadas ya que a cada empleado se le realizan dos mediciones, a saber, el salario con el que entró a la empresa y el salario actual.

Nos debe salir lo siguiente:

Las condiciones para aplicar este test son las siguientes:

Las filas deben ser independientes entre sí (realmente, el requisito es que las diferencias obtenidas al restar un dato con su emparejado sean independientes) Tal y como hemos comentado anteriormente, supondremos que efectivamente es así.
Las diferencias entre los valores pareados deben seguir una distribución normal con media cero. En nuestro caso, la variable construida al restar el salario final del inicial, debe seguir una normal centrada en cero.

Podemos apreciar que la distribución resultante de las diferencias dista mucho de ser una normal (ya que hay muchos puntos que están distantes de la línea de referencia), por lo que aunque el test es bastante robusto a desviaciones de la normalidad, deberíamos aplicar un test no paramétrico.

1.2.2 Test no paramétrico para comparar dos muestras dependientes (medidas repetidas)

dataset: Datos de empleados
test: ¿Puede considerarse que la diferencia entre las medianas del salario inicial y el salario final es cero?

Aplicamos el test de Wilcoxon marcando la casilla correspondiente en el mismo cuadro de diálogo anterior. Debemos obtener lo siguiente:

El test resulta significativo por lo que podemos concluir que, efectivamente, las medianas son distintas.

1.2.3 Anova para k muestras dependientes (test paramétrico de medidas repetidas)

Al igual que en el caso de k muestras independientes, ahora vamos a analizar la situación que se nos presenta cuando tenemos k (>2) muestras dependientes. Por ejemplo, realizamos tres mediciones de una variable sobre cada uno de los sujetos de un experimento y queremos ver si hay diferencias significativas entre las medias (caso paramétrico) o las medianas (caso no paramétrico)

Como el ejemplo de los datos de los empleados no hay variables de este tipo (no se ha registrado el salario en otro momento que no fuese el inicial y el actual), usamos otro conjunto de datos que se puede descargar del siguiente enlace (proviene originalmente de https://maths.shu.ac.uk/mathshelp/SSupport_Practice.html):

Colesterol

Contiene los niveles de colesterol de diferentes individuos. Las mediciones se realizan tres veces, una antes de la aplicación de un tratamiento, otra 4 semanas después y la última 8 semanas después.

1.2.3.1 Test conjunto (omnibus)

dataset: Colesterol
test: ¿Puede considerarse que las medias de las tres mediciones de colesterol son las mismas (k=3)?

En este caso, debemos seleccionar ANOVA/Repeated Measures Anova. En el panel “Repeated Measures Factor” definimos los tres valores o factores que, en este ejemplo, corresponden a los tres instantes de tiempo (“antes”, “después 3 sem”, “después 8 sem”). Cambiamos también la etiqueta “RM Factor 1” por “Tiempo”. A continuación, arrastramos las tres variables al panel “Repeated measures cell”. Como “Dependent variable label” introducimos “Nivel de colesterol”

Nos debe salir lo siguiente:

El test es significativo, por lo que rechazaríamos la hipótesis nula de igualdad en las medias. Así pues, parece que el tratamiento cambia el nivel de colesterol. Como siempre, debemos comprobar las condiciones de aplicabilidad del test. En el caso de medidas repetidas son las siguientes:

Independencia de los datos. Para aplicar el ANOVA de medidas repetidas, al igual que en el caso del Anova de muestras independientes, necesitamos que las filas sean independientes entre sí. Supondremos por tanto que los individuos han sido seleccionados de forma aleatoria. Obviamente, los datos de las columnas no son independientes entre sí (estamos ante un caso de medidas repetidas)
Esfericidad: Las varianzas de las diferencias dos a dos entre todas las columnas de medidas repetidas, han de ser las mismas (el número de comprobaciones obviamente se dispara cuando hay muchas columnas). Para comprobarlo, se aplica el test de Mauchly. Desde Jamovi, hay que seleccionarlo en “Assumption checks / Sphericity tests” y podemos ver que da significativo, por lo que rechazamos que las varianzas de las diferencias sean las mismas. Afortunadamente, para solventar el problema, basta con aplicar una penalización en los grados de libertad. En nuestro caso, desde “Assumption checks / Sphericity corrections” seleccionamos GreenHouse-Geiser y obtenemos lo siguiente:

Como podemos apreciar, aún podemos seguir rechazando la hipótesis nula.
Normalidad de los residuos. Desde el mismo panel de “Assumption Checks” seleccionamos el gráfico QQ y obtenemos:

Vemos que la nube de puntos se distribuye cerca de la línea de referencia, por lo que podemos aceptar que los residuos no se alejan demasiado de una distribución normal.

1.2.3.2 Post hoc

Para realizar las comparaciones múltiples de los tres niveles del factor tiempo, seleccionamos “Post Hoc Tests” y llevamos el factor “Tiempo” al panel de la derecha:

Seleccionamos el método específico de Tukey y el genérico usando Holm como penalización del nivel de significación. Obtenemos lo siguiente:

Podemos apreciar que hay diferencias significativas entre todos los niveles del factor “Tiempo”. Observe que si bien era de esperar entre la primera medición y la segunda y entre la primera y la tercera (hay una diferencia media de 0.5661 y 0.6289 respectivamente) no lo era tanto entre la segunda medición y la tercera (hay una diferencia de tan sólo 0.0628) La razón es que la varianza de los valores de dicha diferencia es muy pequeña, lo que permite afinar el resultado.

1.2.4 Anova para k muestras dependientes (test no paramétrico de medidas repetidas)

Si la hipótesis de normalidad en el anterior test se hubiese violado de una forma muy significativa, tendríamos que aplicar un test no paramétrico como el que vamos a ver en este apartado. Además, en aquellos casos en los que trabajamos con datos ordinales, también debemos aplicar un test de este tipo. Vamos a usar otro dataset ficticio proporcionado por SPSS en el que se registran las preferencias (evaluadas de 1 a 10) por distintos individuos sobre cuatro marcas diferentes (son datos ordinales). Si carga los datos del siguiente enlace:

Preferencias comerciales

obtendrá una vista como la siguiente:

Cada fila corresponde a un individuo y cada columna a una marca. Los valores enteros representan el grado de preferencia de cada individuo a cada marca.

1.2.4.1 Test conjunto (omnibus)

dataset: Preferencias comerciales
test (en el caso de que todas las columnas tengan una distribución similar): ¿Puede considerarse que las medianas de las preferencias es la misma en las cuatro marcas comerciales (k=4)?
test (en el caso de que todas las columnas no tengan una distribución similar): ¿Puede considerarse que ninguna marca domina en sus preferencias sobre el resto (k=4)?

Seleccionamos ANOVA/Repeated measures ANOVA Friedman para aplicar el test de Friedman.

Y podemos observar que el test sale significativo, por lo que podemos rechazar la hipótesis nula y aceptamos por tanto que hay alguna marca que domina sobre el resto. Veamos si se cumplen las hipótesis. En el test de Friedman son las siguientes:

Independencia entre las filas. Tal y como hemos venido asumiendo en las secciones anteriores, supondremos que se cumple.
En el caso de que queramos que el test sea sobre las medianas de cada grupo (es algo más restrictivo), la distribución de cada grupo ha de ser similar. Para comprobarlo, seleccionamos “Exploration / Descriptives”, pasamos las variables al panel de “Variables” y seleccionamos “Histograms / Histogram” y “Histograms / Density”

Podemos considerar que las cuatro variables presentan distribuciones similares por lo que podemos concluir que el test de Friedman es un test sobre las medianas y por tanto podemos afirmar que hay al menos una marca que tiene una mediana distinta al resto.

1.2.4.2 Post hoc

Para realizar el post-hoc, Jamovi sólo ofrece por ahora el test de Durbin: “Pairwise comparisons (Durbin-Conover)” Puede ejecutarlo y apreciar que las medianas que pueden considerarse distintas son “com_1”, “com_2” pero no se puede distinguir entre las dos últimas. Así pues obtendríamos {“com_1”, “com_2”, {“com_3”, “com4”}}

En cualquier caso, Jamovi llama al post hoc de Durbin sin aplicar ninguna penalización (debería aplicar una penalización tipo Holm) y ésto sólo sería correcto si el test omnibus de Durbin (un test alternativo al de Friedman) fuese significativo (consulte la ayuda del paquete PMCMRPlus). Es de esperar que, en futuras versiones, Jamovi afine este apartado. En la parte de R veremos cómo realizar el análisis correcto.

1.3 Test para comparar resultados de varios clasificadores

En este apartado vamos a comparar el comportamiento de varios clasificadores. Puede descargar los datos desde el siguiente enlace:

Classification Algorithms

obtendrá una vista como la siguiente:

El método que vamos a seguir para la comparación es el siguiente: tenemos 7 algoritmos y vamos a medir el porcentaje de acierto en una treintena de conjuntos de datos. Dicho porcentaje ha podido obtenerse, por ejemplo, como una media de los porcentajes de acierto después de una validación cruzada. Queremos ver si podemos considerar que los siete clasificadores tienen una tasa de acierto similar o si hay alguno que no. En este ejemplo, podemos considerar que cada uno de los conjuntos de datos utilizados son los sujetos del estudio, a los que se les ha aplicado una serie de medidas repetidas al evaluar el comportamiento de distintos algoritmos sobre cada uno de ellos. Estamos por tanto en un caso de k muestras dependientes (o de medidas repetidas)

Veamos qué requisitos se cumplen para aplicar un test u otro:

Podemos aceptar que las filas (los datasets) se han elegido de forma aleatoria por lo que son independientes entre sí.
Si analizamos los histogramas de las columnas (Exploration/Descriptives/Histogram) vemos que ninguno se asemeja a una normal. Además, muchos de ellos presentan varios picos lo que supone una desviación importante del requisito de normalidad, por lo que debemos usar un test no paramétrico. Si lo desea, puede ejecutar desde Jamovi un Anova paramétrico y ver que el gráfico QQ de los residuos nos dice que, tal y como era de esperar, éstos no se ajustan a una distribución normal. El test a aplicar, por tanto, es un test no paramétrico. Así pues aplicaremos el test de Friedman (no paramétrico de medidas repetidas)
Una vez que hemos determinado que el test a aplicar es el de Friedman, necesitamos ver si las distribuciones de los grupos (las columnas) son similares. Tal y como hemos apreciado anteriormente, no es así ya que los histogramas eran muy distintos entre sí. Por lo tanto, el test de Friedman es un test de dominancia estocástica y no sobre las medianas.

Si aplicamos el test de Friedman y el post hoc correspondiente, obtenemos los siguientes grupos:

{Alg1, {Alg2, Alg4, Alg5, Alg6}, {Alg3, Alg7}}

Si observamos las tasas de acierto, vemos que Alg1 es un algoritmo muy malo, mientras que Alg3 y Alg7 son mejores que el resto, pero indistiguibles entre sí. En resumen, podemos afirmar que los algoritmos 3 y 7 presentan los mejores resultados, dominando estocásticamente sobre el resto de algoritmos.

Ejecute el Anova paramétrico y observe que en el post hoc sólo obtiene que el algoritmo 1 es distinto al resto. Este es un ejemplo en el que si no se verifican los requisitos del Anova paramétrico no debemos ejecutarlo ya que, en el caso de hacerlo, podríamos obtener (erróneamente) unas conclusiones que no nos permitiesen rechazar la igualdad entre algunos grupos. En igualdad de condiciones, el test de Friedman tiene menos potencia que el Anova paramétrico, pero si no se dan dichas condiciones, simplemente no debemos ejecutar el Anova.

Finalmente, comentemos que en la sección de R, realizaremos un análisis similar pero suprimiremos el Algoritmo 1 ya que presenta sistemáticamente resultados peores y no tiene mucho sentido incluirlo en la comparación. Tenga en cuenta que el test omnibus (el conjunto) obviamente dará significativo y queremos evitar esta situación. Si no incluimos el algoritmo 1 y el test omnibus da significativo, tendremos la garantía de que es al menos uno de los algoritmos decentes el que tiene un comportamiento distinto al resto.

2 Tests de Hipótesis con R

Esta parte es de ampliación, aunque es muy recomendable ya que hay ciertos tests que no se pueden lanzar desde el entorno gráfico de Jamovi. Por ejemplo, por ahora, Jamovi no permite realizar un post-hoc en el que se fije un grupo y por tanto la penalización del nivel de significación sea más suave al disminuir el número de comparaciones múltiples. Esto es de especial importancia cuando se va a comparar el resultado de un clasificador nuevo que hayamos desarrollado con otros clasificadores existentes en la literatura. En esta caso, no queremos comparar los tests ya existentes entre sí sino únicamente nuestro nuevo test (será el control) con el resto.

Muchos de los tests que vamos a ver están directamente accesibles desde R, mientras que otros están en paquetes que habrá que instalar. La lista de paquetes a instalar (y las librerías correspondientes que habrá que cargar en cada sesión) se pueden encontrar en el siguiente enlace:

https://decsai.ugr.es/jccubero/Master/IntrodEstadLibrerias.R

Uno de los paquetes que vamos a usar es rstatx. Proporciona distintos wrappers que permiten usar los tests más conocidos en flujos de ejecución al estilo de tidyverse. También proporciona herramientas para mostrar gráficos de cajas con los p-values obtenidos con un post hoc. Para facilitar su uso, puede utilizar las funciones definidas en el siguiente enlace (© JC Cubero):

https://decsai.ugr.es/jccubero/Master/IntrodEstadFunciones.R

En el anterior paquete puede encontrar las siguientes funciones:

GraficoCajas = function (datos, variable_depend, variable_factor, titulo = "")

Muestra un diagrama de cajas conjunto de una variable dependiente agrupada por los valores de la variable factor.

GraficoCajasTest = function (datos, variable_depend, variable_factor, 
                            es_param, es_mr, test)

Muestra un diagrama de cajas conjunto de una variable dependiente agrupada por los valores de la variable factor, incluyendo los resultados de un test del paquete rstatix. Los parámetros es_param (es paramétrico) y es_mr (es de medidas repetidas) son lógicos y se usan para construir el título del gráfico.

GraficoCajasPost = function(datos, variable_depend, variable_factor, 
                            es_param, es_mr, test, 
                            control = "", estrellas = FALSE)

Realiza lo mismo que la anterior función para un test post hoc, por lo que añade los p-values de las comparaciones múltiples. Si hay un grupo de control se pasa en el parámetro control. El parámetro estrellas indica si se quieren imprimir asteriscos en vez de los valores numéricos de los p-values

Para ejecutar los distintos tests, puede descargar los siguientes conjuntos de datos:

Datos de empleados (RData) contiene el conjunto de datos con nombre DatosEmpleados
Cholesterol (RData) contiene el conjunto de datos con nombre Colesterol

Al igual que en la la sección tests con Jamovi suponemos que las filas son independientes entre sí y que las muestras son independientes dentro de cada grupo (también entre los distintos grupos en el caso de muestras independientes).

2.1 Tests para muestras independientes

Vamos a cargar el dataset DatosEmpleados y establecemos los nombres de las variables con las que vamos a trabajar. También convertimos a factor la columna de la variable de agrupación ya que algunos tests como aov así lo requieren. Lo hacemos con la función as.factor:

load("DatosEmpleados.RData")
datos = DatosEmpleados
v.depend = "salario"
v.factor.2.valores = "sexo"
indice.col.factor.2valores = which(colnames(datos) == v.factor.2.valores)
datos[,indice.col.factor.2valores] = as.factor(as.data.frame(datos)[,indice.col.factor.2valores])

head(datos)
##   Id   sexo   fechnac educ         catlab salario salini tiempemp expprev minoría
## 1  1 Hombre  2/3/1952   15      Directivo   57000  27000       98     144      No
## 2  2 Hombre 5/23/1958   16 Administrativo   40200  18750       98      36      No
## 3  3  Mujer 7/26/1929   12 Administrativo   21450  12000       98     381      No
## 4  4  Mujer 4/15/1947    8 Administrativo   21900  13200       98     190      No
## 5  5 Hombre  2/9/1955   15 Administrativo   45000  21000       98     138      No
## 6  6 Hombre 8/22/1958   15 Administrativo   32100  13500       98      67      No

Algunos tests usan fórmulas para especificar la variable dependiente y la de agrupación. Por ejemplo, si queremos realizar un Anova del salario sobre el sexo, usaríamos la fórmula siguiente:

salario ~ sexo

Para no tener que repetir dicha expresión, la guardamos en una variable formula aplicando la función as.formula:

formula.t.test = as.formula(paste0(v.depend, "~", v.factor.2.valores))  # 2 muestras independientes

formula.t.test
## salario ~ sexo

2.1.1 Prueba para dos muestras independientes (test paramétrico)

Trabajamos sobre DatosEmpleados con las variables v.depend (salario) y v.factor.2.valores (sexo). En primer lugar, construimos un gráfico con los diagramas de cajas:

GraficoCajas(datos, v.depend, v.factor.2.valores, "Boxplot")

El gráfico nos muestra que hay mucha menos varianza en el grupo de las mujeres que en el de hombres. Para confirmarlo podemos llamar a la función levene_test. El test de Levene sale, claramente, significativo, por lo que rechazamos que las varianzas sean iguales:

test_homoc_levene = levene_test(formula.t.test, data = datos)
test_homoc_levene
## # A tibble: 1 × 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1     1   472      51.6 2.67e-12

Como las varianzas son distintas, lanzamos el test de Welch llamando a la correspondiente función de rstatix y mostramos el resultado con la función GraficoCajasTest.

test_welch = welch_anova_test(formula = formula.t.test, data = datos)

test_welch contiene información sobre el valor de estadístico del test, los grados de libertad, el p-value, etc. Para mostrar este último, basta acceder al valor p:

test_welch$p
## [1] 8.3e-27

O directamente lo vemos en la gráfica correspondiente (en vez de mostrar el valor exacto, simplemente nos dice que es menor que 0.0001):

GraficoCajasTest(datos, v.depend, v.factor.2.valores, 
                 es_param = TRUE, es_mr = FALSE, test_welch)

Aunque el test de Welch es robusto a la desviaciones de la normalidad de los residuos, debemos analizarlos para tener una garantía estadística. Para extraer los residuos de un test realizado con rstatix, debemos utilizar la función anova_test que veremos posteriormente (es la generalización de la función t_test cuando hay más de dos grupos) Definimos nuestra propia función Residuos para realizar dicha tarea:

Residuos = function (data, var.depend, var.factor){
  formula.anova = as.formula(paste0(var.depend, "~", var.factor))  # k>2 muestras independientes
  test_anova = anova_test(formula = formula.anova, data = data)
  modelo_test = attr(test_anova, "args")$model
  residuals(modelo_test)
}

test_anova_2muestras = anova_test(formula.t.test, data = datos)
residuos_test_anova_2muestras = Residuos(datos, v.depend, v.factor.2.valores)

Para ver si siguen una distribución normal, podemos lanzar el test de Shapiro aunque normalmente es suficiente con que el gráfico QQ no muestre una desviación acusada (usamos las funciones shapiro.test y ggqqplot)

test_normalidad_residuos = shapiro.test(x = residuos_test_anova_2muestras)
test_normalidad_residuos
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_test_anova_2muestras
## W = 0.83621, p-value < 2.2e-16
ggqqplot(residuos_test_anova_2muestras)

Podemos apreciar (ya lo habíamos visto con Jamovi) que no se cumple la hipótesis de normalidad de los residuos, por lo que optamos por aplicar un test no paramétrico. Lo hacemos a continuación.

2.1.2 Test no paramétrico para comparar dos muestras independientes

Para aplicar el test de Mann-Whitney llamamos a la función wilcox_test de rstatix y construimos la gráfica llamando a la función GraficoCajasTest:

test_mann_whitney = wilcox_test(data = datos, formula = salario ~ sexo)
GraficoCajasTest(datos, v.depend, v.factor.2.valores, 
                 es_param = TRUE, es_mr = FALSE, test_mann_whitney)

Nota: El valor del estadístico \(W\) obtenido por la función wilcox_test (\(46111.5\)) es el mayor entre \(U_1\) y \(U_2\) (consultar transparencias del curso) mientras que Jamovi muestra el minimo (\(9617\))

Recordemos que el test de Mann-Whitney es un test de dominancia estocástica. Si los grupos presentan distribuciones similares, también será un test sobre las medianas, pero ya habíamos visto que las varianzas del salario eran distintas en el grupo de hombres y en el grupo de mujeres, por lo que las distribuciones no son similares. Por tanto, no podemos afirmar que sea un test sobre las medianas y concluimos que el salario de los hombres es una variable tal que si elegimos al azar un hombre y una mujer es más probable que el salario del primero sea mayor que el de la mujer.

2.1.3 Anova para k muestras independientes (test paramétrico)

2.1.3.1 Test conjunto (omnibus)

Trabajamos sobre DatosEmpleados con las variables salario y catlab. En primer lugar, establecemos los nombres de las variables y de la fórmula:

v.factor = "catlab"
indice.col.factor = which(colnames(datos) == v.factor)
datos[,indice.col.factor] = as.factor(as.data.frame(datos)[,indice.col.factor]) 
formula.anova = as.formula(paste0(v.depend, "~", v.factor))  # k>2 muestras independientes

Para comprobar la normalidad de los residuos, llamamos a la función Residuos que habíamos definido previamente:

residuos.anova = Residuos (data = datos, v.depend, v.factor)
ggqqplot(residuos.anova)

Construimos un gráfico con los diagramas de cajas y aplicamos el test de Levene:

test_homoc_levene = levene_test(formula.anova, data = datos)
test_homoc_levene
## # A tibble: 1 × 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1     2   471      51.2 7.65e-21

GraficoCajas(datos, v.depend, v.factor, "Boxplot")

El test de Levene confirma lo que vemos en el diagrama de cajas y es que las varianzas son distintas en cada grupo. Aplicamos por tanto el test de Welch y nos ha de salir lo siguiente:

# COMPLETAR

El test es significativo y por tanto podemos asegurar que hay al menos una categoría laboral con un salario distinto al resto. Procedemos al post-hoc:

2.1.3.2 Post hoc

Post hoc específico

Como los grupos según la categoría laboral tienen varianzas distintaS, Lanzamos las comparaciones múltiples con el post-hoc específico de Games-Howell llamando a la función games_howell_test de rstatix con los parámetros habituales de la fórmula y el conjunto de datos. Llamamos también a la función GraficoCajasPost y vemos que el test rechaza la igualdad entre las medias de todos los grupos:

post_games_howell <- games_howell_test(datos, formula.anova)     
GraficoCajasPost(datos, v.depend, v.factor, 
                 es_param = TRUE , es_mr = FALSE, post_games_howell)

Tal y como vimos en la sección de Jamovi, si hubiésemos utilizado el test de Tukey HSD como post hoc, no podríamos haber rechazado la igualdad entre Administrativos y Seguridad ya que es test de Tukey sólo ha de aplicarse cuando todos los grupos tienen la misma varianza y en este caso no se cumple. Para comprobarlo, lanzamos el test con la función tukey_hsd de rstatix:

post_tukey <- tukey_hsd(datos, formula.anova)     
GraficoCajasPost(datos, v.depend, v.factor, 
                 es_param = TRUE , es_mr = FALSE, post_tukey)

Post hoc genérico

Vamos a ejecutar ahora las comparaciones múltiples usando para cada comparación un test generico y aplicando una penalización. En este caso, usaríamos el test de Welch para realizar cada una de las comparaciones entre las medias de dos grupos. Recordemos que Jamovi (apartado Anova para k muestras independientes (test paramétrico)) usaba en el test genérico una varianza agregada pero no era lo recomendable al tener los grupos varianzas diferentes. Por tanto, nuestro propósito es lanzar las comparaciones múltiples, usando como test base para cada comparación el test de Welch (no usando la varianza agregada) y penalizando con Holm. Para ello, llamamos a la función pairwise_t_test de rstatix (el parámetro pool.sd es el de la varianza agregada y es TRUE por defecto por lo que lo tenemos que cambiar)

post_pw_t_test <- pairwise_t_test (formula.anova, data = datos,
                                   pool.sd = FALSE)
GraficoCajasPost(datos, v.depend, v.factor , 
                 es_param = TRUE , es_mr = FALSE, post_pw_t_test)

Si queremos obtener el mismo resultado que Jamovi (algo no recomendado en este caso por lo que hemos comentado anteriormente) habríamos ejecutado lo siguiente:

post_pw_t_test_JAMOVI <- pairwise_t_test (formula.anova, data = datos,
                                   pool.sd = TRUE)
post_pw_t_test_JAMOVI
## # A tibble: 3 × 9
##   .y.     group1         group2       n1    n2         p p.signif     p.adj p.adj.signif
## * <chr>   <chr>          <chr>     <int> <int>     <dbl> <chr>        <dbl> <chr>       
## 1 salario Administrativo Directivo   363    84 9.61e-109 ****     2.88e-108 ****        
## 2 salario Administrativo Seguridad   363    27 1.26e-  1 ns       1.26e-  1 ns          
## 3 salario Directivo      Seguridad    84    27 1.27e- 40 ****     2.54e- 40 ****

Nota. También podríamos haber llamado a la función t_test de rstatix (la que usábamos para comparar dos muestras) con el parámetro var.equal = FALSE (utilizará la varianza no agregada) y p.adjust.method="holm" ya que si le pasamos a t_test como parámetro una fórmula que tiene un factor con más de dos niveles, automáticamente realiza las comparaciones múltiples llamando a pairwise_t_test

Post hoc con control

Si deseamos fijar un grupo de control y así disminuir la penalización para corregir el error FWER, podemos lanzar el test específico de Dunnett disponible en el paquete DescTools. También podemos usar el mismo test genérico usando la función pairwise_t_test pasándolo el grupo de control en el parámetro ref.group

2.1.4 Anova para k muestras independientes (test no paramétrico)

2.1.4.1 Test conjunto (omnibus)

Recordemos que al no cumplirse la hipótesis de normalidad de los residuos, es conveniente aplicar un test no paramétrico. Para ejecutar el test de Kruskal-Wallis llamamos a la función kruskal_test de rstatix en la forma habitual con la fórmula y los datos. Llamamos también a la función GraficoCajasTest. Nos debe salir lo siguiente:

# COMPLETAR

Recordemos que el test de Kruskal-Wallis es un test de dominancia estocástica. Si los grupos hubiesen presentado distribuciones similares, también sería un test sobre las medianas. Sin embargo, hemos visto que las varianzas son distintas en cada grupo por lo que no tienen distribuciones similares. En cualquier caso, si queremos ver las distribuciones de cada grupo, podemos usar la función ggdensity de la siguiente forma:

ggdensity(as.data.frame(datos), v.depend, facet.by = v.factor)

2.1.4.2 Post hoc

Post hoc específico

rstatix proporciona el test específico de Dunn en vez del test de Dwass-Steel-Critchlow-Fligner que veíamos con Jamovi. Se ejecuta en la forma habitual:

post_kw_dunn = dunn_test(formula.anova, data = datos)
GraficoCajasPost(datos, v.depend, v.factor,  
                 es_param = FALSE , es_mr = FALSE, post_kw_dunn)

Un test específico aún mejor que el de Dunn es el test de Conover incluido en el paquete PMCMRPlus a través de la función kwAllPairsConoverTest. En este test, hay que especificar el tipo de penalización a realizar. Usamos la de Holm:

post_kw_conover = kwAllPairsConoverTest(formula.anova, data = datos,
                                        p.adjust.method = "holm")
post_kw_conover
##           Administrativo Directivo
## Directivo < 2e-16        -        
## Seguridad 1.9e-05        3.1e-10

Podemos apreciar que el test de Conover es capaz de rechazar la igualdad entre Administrativos y Seguridad con más fuerza que el de Dunn.

Nota. Si desea lanzar manualmente el test de Dwass-Steel-Critchlow-Fligner, puede recurrir a la función dscfAllPairsTest en el paquete PMCMRPlus. También puede ejecutar el test de Nemenyi a través de la función kwAllPairsNemenyiTest del mismo paquete.

Nota. EL paquete PMCMRPlus contiene muchos tests para realizar el post hoc no paramétrico. Sin embargo, los tests no son compatibles con rstatix, por lo que no podremos mostrar los resultados con la función GraficoCajasPost.

Post hoc genérico

Para ejecutar las comparaciones múltiples usando un test genérico con una penalización del nivel de significación usamos la misma función wilcox_test que usábamos para ejecutar el test de Mann-Whitney en el caso no paramétrico de 2 muestras independientes. Si la variable factor de la fórmula pasada como parámetro contiene más de un grupo, automáticamente se realizan las comparaciones múltiples:

post_kw_wilc = wilcox_test(formula.anova, data = datos)
GraficoCajasPost(datos, v.depend, v.factor, 
                 es_param = FALSE, es_mr = FALSE, post_kw_wilc)

Concluimos por tanto que hay diferencias significativas en los tres grupos. Formalmente, para poder concluir cuál es mayor que la otra, deberíamos haber aplicado tests a una cola (one-sided test) Para ello habría que seleccionar la opción alternative = "greater" en la función wilcox-test. En cualquier caso, es claro que los directivos dominan de forma estocástica sobre los otros dos grupos. No está tan claro en el caso de los de Seguridad y Administrativos, por lo que si quisiéramos inferir conclusiones con garantía estadística, habría que lanzar los tests a una cola (si lo hace comprobará que, efectivamente, el salario del grupo de seguridad es mayor -según el criterio de dominancia estocástica- que los administrativos)

Post hoc con control

Si queremos fijar un grupo de control, podemos usar el test específico de Conover (usando la función kwManyOneConoverTest de PMCMRPlus) o bien el test genérico de Wilcoxon usando la función anterior wilcox_test y el parámetro ref.group.

2.2 Tests para muestras dependientes o de medidas repetidas

Tal y como vimos en el apartado de Jamovi, asumimos que las filas dentro de cada grupo son independientes entre sí. Obviamente, las filas de grupos distintos no son independientes entre sí.

Usaremos el sufijo mr (medidas repetidas) en los nombres de las variables y objetos que creemos.

2.2.1 Prueba T (t-Student) para dos muestras pareadas (medidas repetidas)

test: ¿Puede considerarse que la diferencia entre las medias del salario actual y el salario inicial es cero?

Para lanzar el test T para muestras pareadas usando directamente el conjunto de datos original, podemos usar la función t.test (observe que no es la función t_test de rstatix) Le pasamos como parámetro los dos vectores correspondientes a las dos muestras pareadas (salario inicial y salario final) y, muy importante, el parámetro paired = TRUE (muestras dependientes):

col.salario.inicial = datos$salario
col.salario.final = datos$salini
t.test.mr = t.test(col.salario.inicial, col.salario.final, 
                   paired = TRUE)
t.test.mr
## 
##  Paired t-test
## 
## data:  col.salario.inicial and col.salario.final
## t = 35.036, df = 473, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  16427.41 18379.56
## sample estimates:
## mean difference 
##        17403.48

En cualquier caso, la mayor parte de los tests sobre medidas repetidas los lanzaremos de otra forma. Necesitamos convertir los datos desde el formato wide al formato long. En el ejemplo de los datos de empleados, queremos construir un nuevo conjunto de datos con dos variables nuevas. Por una parte tendremos la variable nominal Tiempo que puede verse como un factor intra sujetos. Dicha variable tendrá dos valores de tipo texto: "salini" (salario inicial) y "salario" (salario final). Observe que ambos valores son nominales (cadenas de texto) que coindicen con los nombres de las columnas originales. Por otra parte, tendremos una nueva variable Salario (se ha usado una mayúscula al inicio para distinguirla de la columna original salario). Cada fila del conjunto original se desdoblará en dos: una fila contendrá el valor (texto) "salini" en la columna Tiempo y el valor (dato de tipo real) del salario inicial en la columna Salario. La otra fila contendrá el valor de texto "salario" en la columna Tiempo y el valor del salario final (dato de tipo real) en la columna Salario. Este mismo proceso se generaliza con k grupos.

Para realizar esta conversión, hacemos uso de la función pivot_longer del paquete tidyverse. Almacenamos en la variable v.id.mr el nombre de la columna del conjunto de datos que identifica a cada fila de forma unívoca (en nuestro caso Id). En el vector columnas.mr (columnas.medidas.repetidas) guardamos los nombres de las columnas que representan las medidas repetidas. En la variable v.depend.mr guardamos el nombre de la nueva variable dependiente que se va a crear. En la variable v.factor.mr guardamos el nombre de la nueva variable factor que vamos a crear, en nuestro caso Tiempo. El resultado lo guardamos en el dataset datos.mr.long:

datos.mr.wide = DatosEmpleados
v.id.mr = "Id"
columnas.mr = c("salario", "salini")
v.depend.mr = "Salario"
v.factor.mr = "Tiempo"
datos.mr.long = pivot_longer(datos.mr.wide, 
                             cols = any_of(columnas.mr), 
                             names_to = v.factor.mr,
                             values_to = v.depend.mr)

datos.mr.wide[c(1:3),]
##   Id   sexo   fechnac educ         catlab salario salini tiempemp expprev minoría
## 1  1 Hombre  2/3/1952   15      Directivo   57000  27000       98     144      No
## 2  2 Hombre 5/23/1958   16 Administrativo   40200  18750       98      36      No
## 3  3  Mujer 7/26/1929   12 Administrativo   21450  12000       98     381      No
datos.mr.long[c(1:6),]
## # A tibble: 6 × 10
##      Id sexo   fechnac    educ catlab         tiempemp expprev minoría Tiempo  Salario
##   <int> <fct>  <chr>     <int> <fct>             <int>   <int> <fct>   <chr>     <int>
## 1     1 Hombre 2/3/1952     15 Directivo            98     144 No      salario   57000
## 2     1 Hombre 2/3/1952     15 Directivo            98     144 No      salini    27000
## 3     2 Hombre 5/23/1958    16 Administrativo       98      36 No      salario   40200
## 4     2 Hombre 5/23/1958    16 Administrativo       98      36 No      salini    18750
## 5     3 Mujer  7/26/1929    12 Administrativo       98     381 No      salario   21450
## 6     3 Mujer  7/26/1929    12 Administrativo       98     381 No      salini    12000

Convertimos también al tipo factor la variable nueva factor que hemos creado, así como la columna del identificador. Lo hacemos con los nombres genéricos que hemos introducido anteriormente:

indice.col.id.mr = which(colnames(datos.mr.long) == v.id.mr)
indice.col.factor.mr = which(colnames(datos.mr.long) == v.factor.mr)
datos.mr.long[,indice.col.id.mr] = as.factor(as.data.frame(datos.mr.long)[,indice.col.id.mr])
datos.mr.long[,indice.col.factor.mr] =as.factor(as.data.frame(datos.mr.long)[,indice.col.factor.mr])

Ya tenemos preparado nuestro conjunto de datos para aplicar las funciones de rstatix (y de muchos otros paquetes de R). Si vamos a usar la función t_test (observe que es t_test y no t.test) tenemos que usar la siguiente fórmula:

formula.mr.sin.error = Salario ~ Tiempo

o usando el formato genérico con las variables anteriores:

formula.mr.sin.error = as.formula(paste0(v.depend.mr, "~", v.factor.mr))
formula.mr.sin.error
## Salario ~ Tiempo

Llamamos a la función t_test con el parámetro paired = TRUE

test_t_mr <- t_test(formula = formula.mr.sin.error, 
                    data = datos.mr.long, paired = TRUE) 
GraficoCajasTest(datos.mr.long, v.depend.mr, v.factor.mr,  
                 es_param = TRUE, es_mr = TRUE, test_t_mr)

Para comprobar la hipótesis de normalidad basta con calcular la columna formada por las diferencias de las dos muestras pareadas (ambas tienen el mismo tamaño) y aplicar el test de Shapiro o ver el gráfico QQ. Podemos usar las columnas que habíamos construido (col.salario.inicial y col.salario.final) o bien creamos unas nuevas con los nombres genéricos que hemos introducido anteriormente:

col.mr.1 = datos.mr.wide[,columnas.mr[1]]
col.mr.2 = datos.mr.wide[,columnas.mr[2]]
shapiro.test(x = col.mr.1 - col.mr.2)
## 
##  Shapiro-Wilk normality test
## 
## data:  col.mr.1 - col.mr.2
## W = 0.78168, p-value < 2.2e-16
ggqqplot(col.mr.1 - col.mr.2)

Podemos apreciar que hay una desviación bastante notable de la normalidad por lo que será preferible aplicar un test no paramétrico.

Nota. También podríamos haber analizado los residuos: esto se verá posteriormente, cuando veamos el ANOVA general para medidas repetidas.

2.2.2 Test no paramétrico para comparar dos muestras dependientes (medidas repetidas)

test: ¿Puede considerarse que la diferencia entre las medianas del salario inicial y el salario final es cero?

En el caso no paramétrico, debemos lanzar el test de Wilcoxon (Wilcoxon signed-rank test). Para ello, llamamos a la función wilcox_test de rstatix. Le pasamos la misma fórmula que habíamos usado en el test t de 2 muestras con medidas repetidas e incluimos el parámetro paired = TRUE:

formula.mr.sin.error
## Salario ~ Tiempo
test_wilcoxon_mr = wilcox_test(formula.mr.sin.error, data = datos.mr.long, 
                               paired = TRUE)
GraficoCajasTest(datos.mr.long, v.depend.mr, v.factor.mr,
                 es_param = FALSE, es_mr = TRUE ,
                 test_wilcoxon_mr)

Para poder asegurar que el test de Wilcoxon es fiable necesitamos comprobar el requisito de que las diferencias se distribuyen de forma simétrica con respecto a la mediana.

diferencias = col.mr.1 - col.mr.2
ggdensity(data = diferencias)

median(diferencias)
## [1] 14250

Vemos que no es una distribución simétrica, por lo que procedemos a aplicar el sign-test aplicando la función sign_test de rstatix. El sign-test sólo tiene en cuenta si un valor es mayor que otro, pero no cuánto más, por lo que pierde información y tiene menos potencia que el test de Wilcoxon.

test_sign_mr = sign_test(formula.mr.sin.error, data = datos.mr.long)

GraficoCajasTest(datos.mr.long, v.depend.mr, v.factor.mr,
                 es_param = FALSE, es_mr = TRUE ,
                 test_sign_mr)

Por lo tanto, a pesar de haber tenido que aplicar el sign-test que es menos potente que el de Wilcoxon, podemos concluir que las medianas son distintas y por tanto el salario mediano inicial es menor que el salario mediano final. Formalmente, para poder concluir cuál es mayor que la otra, deberíamos haber aplicado un test a una cola (one-sided test) Para ello habría que seleccionar la opción alternative = "greater" en la función wilcox_test (o en la función sign_test).

2.2.3 Anova para k muestras dependientes (test paramétrico de medidas repetidas)

Trabajamos con el conjunto de datos Cholesterol.RData. Los datos vienen en formato wide por lo que tendrá que transformarlos en long. Tenga en cuenta que el nombre del dataset es Cholesterol, el vector de columnas de medidas repetidas es c("Before", "After4weeks", "After8weeks"), el nombre de la nueva variable dependiente puede ser NivelColesterol, el nombre de la nueva variable factor puede ser Tiempo y el nombre de la columna de identificador es ID.

load("Cholesterol.RData")
datos.mr.wide = Cholesterol
nombres.columnas = colnames(datos.mr.wide)  
v.id.mr = nombres.columnas[1]  # El identific. es la primera columna
v.depend.mr = "NivelColesterol"
v.factor.mr = "Tiempo"
columnas.mr = nombres.columnas[-1]  # El resto de columnas son las MR

datos.mr.long = pivot_longer(datos.mr.wide, 
                             cols = any_of(columnas.mr), 
                             names_to = v.factor.mr,
                             values_to = v.depend.mr)

indice.col.id.mr = which(colnames(datos.mr.long) == v.id.mr)
indice.col.factor.mr = which(colnames(datos.mr.long) == v.factor.mr)
datos.mr.long[,indice.col.id.mr] = as.factor(as.data.frame(datos.mr.long)[,indice.col.id.mr])
datos.mr.long[,indice.col.factor.mr] =as.factor(as.data.frame(datos.mr.long)[,indice.col.factor.mr])

formula.mr.sin.error = as.formula(paste0(v.depend.mr, "~", v.factor.mr))

Obtenemos los siguientes datos (mostramos como quedarían las dos primeras filas:

datos.mr.wide[c(1:2),]
##   ID Before After4weeks After8weeks
## 1  1   6.42        5.83        5.75
## 2  2   6.76        6.20        6.13
datos.mr.long[c(1:6),]
## # A tibble: 6 × 3
##   ID    Tiempo      NivelColesterol
##   <fct> <fct>                 <dbl>
## 1 1     Before                 6.42
## 2 1     After4weeks            5.83
## 3 1     After8weeks            5.75
## 4 2     Before                 6.76
## 5 2     After4weeks            6.2 
## 6 2     After8weeks            6.13

2.2.3.1 Test conjunto (omnibus)

Para contemplar el caso de k grupos de medidas repetidas, tenemos que tabajar con la función anova_test y con la siguiente fórmula

formula.anova.mr.error = NivelColesterol ~ Tiempo + Error(ID / Tiempo)

Para poder utilizar los nombres genéricos que hemos introducido en los apartados anteriores, usamos la siguiente expresión general:

formula.anova.mr.error = as.formula(paste0(v.depend.mr, "~", 
                                    v.factor.mr, "+ Error(" , 
                                    v.id.mr , "/" , v.factor.mr , ")"))
formula.anova.mr.error
## NivelColesterol ~ Tiempo + Error(ID/Tiempo)

Ahora ya estamos en condiciones de aplicar la función anova_test de rstatix. Dicha función analiza el formato de la fórmula y al encontrar el término Error aplica automáticamente un ANOVA de medidas repetidas:

test_anova_mr <- anova_test(formula = formula.anova.mr.error, 
                            data = datos.mr.long) 
GraficoCajasTest(datos.mr.long, v.depend.mr, v.factor.mr,  
                 es_param = TRUE, es_mr = TRUE, test_anova_mr)

Debemos señalar que la función anova_test automáticamente aplica el test de Mauchly para comprobar la hipótesis de esfericidad y en el caso de que no se cumpla, aplica también automáticamente la corrección de Greenhouse-Geisser

Como ejercicio, escriba la fórmula correspondiente al ejemplo de los datos de los empleados y lance la función anova_test. El resultado es equivalente al que obtuvimos en el apartado test para dos muestras dependientes con la función t_test (puede comprobar que el valor de F obtenido en el ANOVA es el cuadrado de la t)

Para ver la normalidad de los residuos, definimos la siguiente función que extrae los residuos del modelo. Al igual que vimos en el apartado Prueba para dos muestras independientes (test paramétrico), usamos una función similar (hemos definido ResiduosMR usando la función anova_test de rstatix)

ResiduosMR = function(data.long, var.depend.mr, var.factor.mr, var.id.mr, columnas.mr){
  formula.anova.mr.error = as.formula(paste0(var.depend.mr, "~", 
                                    var.factor.mr, "+ Error(" , 
                                    var.id.mr , "/" , var.factor.mr , ")"))
  test_anova_mr = anova_test(formula.anova.mr.error, data = data.long)
  residuos_anova_test_mr = residuals(attr(test_anova_mr, "args")$model)
  # Me da los residuos por cada columna, por lo que hay que agregarlos todos:
  residuos.juntos.mr = pivot_longer(as.data.frame(residuos_anova_test_mr), 
                                    cols = any_of(columnas.mr))
  residuos.juntos.mr  = residuos.juntos.mr$value
}

Llamamos a la función y mostramos los residuos:

formula.anova.mr.error
## NivelColesterol ~ Tiempo + Error(ID/Tiempo)
residuos.mr = ResiduosMR(datos.mr.long, v.depend.mr, v.factor.mr,
                         v.id.mr, columnas.mr)
ggqqplot(residuos.mr)

Podemos apreciar que los residuos se ajustan a una normal, por lo que se cumplen las hipótesis del test paramétrico (realmente no se cumplía el requisito de esfericidad pero el test lo solventa automáticamente aplicando la corrección sobre los grados de libertad de Greenhouse-Geisser) Por tanto podemos concluir que existe al menos un grupo con una media distinta de los otros grupos: en nuestro ejemplo, el nivel de colesterol al inicio del estudio es mayor que el nivel medido después de 4 y 8 semanas de tratamiento. A continuación, procedemos a ejecutar el post-hoc para corroborar la anterior afirmación y para comprobar si hay diferencias significativas entre el nivel de colesterol medido a las 4 y a las 8 semanas.

2.2.3.2 Post hoc

Post hoc específico

Para lanzar un test específico de comparaciones múltiples como el de Tukey, tendríamos que recurrir a otros paquetes y ajustar otros modelos que se escapan del nivel de este guion. Para más información, puede consultar:

https://stackoverflow.com/questions/74839746/anova-with-repeated-measures-and-tukeyhsd-post-hoc-test-in-r

Post hoc genérico

Para lanzar las comparaciones múltiples con un test genérico y una penalización del nivel de significación (por ejemplo, Holm) usamos la función pairwise_t_test de rstatix (es un wrapper de pairwise.t.test). Es muy importante que observe que se lanza con la misma fórmula que habíamos usado en el caso de 2 muestras dependientes y no la fórmula que contenía el término del error. También incluimos el parámetro paired = TRUE:

formula.mr.sin.error = as.formula(paste0(v.depend.mr, "~", v.factor.mr))
formula.mr.sin.error
## NivelColesterol ~ Tiempo

post_pw_t_test_mr = pairwise_t_test(formula.mr.sin.error, 
                                    data = datos.mr.long,
                                    paired = TRUE, 
                                    p.adjust.method = "holm")
GraficoCajasPost(datos.mr.long, v.depend.mr, v.factor.mr, 
                 es_param = TRUE , es_mr = TRUE,
                 post_pw_t_test_mr)

El resultado del test nos dice que también hay diferencias significativas entre el primer y segundo grupo (medición después de 4 y 8 semanas) Si lanzamos el mismo test planteando como hipótesis alternativa “la media de un grupo es mayor o igual que la del otro grupo” (opción alternative = "greater" de la función pairwise_t_test), obtendríamos una significación de 0.002 entre los dos primeros grupos (sigue siendo significativa) En definitiva, podemos concluir que gracias al tratamiento se aprecia una bajada del nivel de colesterol muy acusada en las 4 primeras semanas y otra bajada, mucho menor pero significativa, entre la cuarta y octava semana.

Post hoc con control

Si queremos fijar un grupo de control, podemos llamar a la misma función anterior pairwise_t_test, indicándolo en el parámetro ref.group

2.2.4 Anova para k muestras dependientes (test no paramétrico de medidas repetidas)

Aplicaremos este tipo de tests cuando no se de el requisito de normalidad del caso paramétrico o cuando los datos sean ordinales, como es el caso del ejemplo sobre preferencias comerciales que vimos en la sección de Jamovi: Friedman. Vamos a trabajar con este conjunto de datos. La versión en R la puede descargar del siguiente enlace:

Preferencias comerciales (RData)

Establecemos los nombres de las variables y realizamos la transformación al formato long:

load("CommercialRatings.RData")
datos.mr.wide = CommercialRatings
nombres.columnas = colnames(datos.mr.wide)  
v.id.mr = nombres.columnas[1]  # El identific. es la primera columna
v.depend.mr = "Valoracion"
v.factor.mr = "Marca"
columnas.mr = nombres.columnas[-1]  # El resto de columnas son las MR

datos.mr.long = pivot_longer(datos.mr.wide, 
                             cols = any_of(columnas.mr), 
                             names_to = v.factor.mr,
                             values_to = v.depend.mr)

indice.col.id.mr = which(colnames(datos.mr.long) == v.id.mr)
indice.col.factor.mr = which(colnames(datos.mr.long) == v.factor.mr)
datos.mr.long[,indice.col.id.mr] = as.factor(as.data.frame(datos.mr.long)[,indice.col.id.mr])
datos.mr.long[,indice.col.factor.mr] =as.factor(as.data.frame(datos.mr.long)[,indice.col.factor.mr])

datos.mr.wide[c(1:2),]
##   id com_1 com_2 com_3 com_4
## 1  1     2     6     7     6
## 2  2     1     7     7     6

datos.mr.long[c(1:8),]
## # A tibble: 8 × 3
##   id    Marca Valoracion
##   <fct> <fct>      <int>
## 1 1     com_1          2
## 2 1     com_2          6
## 3 1     com_3          7
## 4 1     com_4          6
## 5 2     com_1          1
## 6 2     com_2          7
## 7 2     com_3          7
## 8 2     com_4          6

2.2.4.1 Test conjunto (omnibus)

El test aplicable en este caso es el de Friedman. La fórmula que debemos usar ahora es de la forma:

var.dependiente ~ factor | identificador
## var.dependiente ~ factor | identificador

En nuestro caso sería:

Valoracion ~ Marca | id
## Valoracion ~ Marca | id

o en general, usando los nombres de las variables que habíamos creado:

formula.anova.friedman_test = as.formula(paste0(v.depend.mr, "~", 
                                                v.factor.mr , "|" , 
                                                v.id.mr))

Basta llamar a la funcion friedman_test de rstatix:

formula.anova.friedman_test
## Valoracion ~ Marca | id

test_friedman = friedman_test (formula.anova.friedman_test, 
                               data = datos.mr.long)
GraficoCajasTest(datos.mr.long, v.depend.mr, v.factor.mr,
                 es_param = FALSE, es_mr = TRUE ,
                 test_friedman)

Finalmente, para comprobar si todos los grupos tienen distribuciones similares, lanzamos la función ggdensity:

Al tener distribuciones similares, podemos afirmar que el test de Friedman es un test sobre las medianas (en caso contrario sería de dominancia estocástica) Así pues, podemos concluir que hay al menos una marca comercial que tiene una preferencia distinta a las otras marcas.

2.2.4.2 Post hoc

Post hoc específico

Jamovi lanza el test de Durbin, pero como ya vimos, no aplica ninguna penalización interna. Lo correcto sería lanzarlo con una corrección como Holm. Para lanzar un test específico de comparaciones múltiples para Friedman, podemos recurrir a los tests incluidos en el paquete PMCMRPlus (Nemenyi, Durbin, etc.) Muchos de los tests de PMCMRPlus requieren que se indiquen los siguientes parámetros:

y es un vector con los datos de la variable dependiente
groups es un vector con los datos del factor
blocks es un vector con los datos del identificador

Los vectores anteriores se refieren al conjunto de datos en el formato long. Construimos por tanto los siguientes vectores:

col.v.depend.mr = as.data.frame(datos.mr.long)[,v.depend.mr]
col.v.factor.mr = as.data.frame(datos.mr.long)[,v.factor.mr]
col.id.factor.mr = as.data.frame(datos.mr.long)[,v.id.mr]

head(datos.mr.long)
## # A tibble: 6 × 3
##   id    Marca Valoracion
##   <fct> <fct>      <int>
## 1 1     com_1          2
## 2 1     com_2          6
## 3 1     com_3          7
## 4 1     com_4          6
## 5 2     com_1          1
## 6 2     com_2          7

a_mostrar = c(1:6)
col.v.depend.mr[a_mostrar]
## [1] 2 6 7 6 1 7
col.v.factor.mr[a_mostrar]
## [1] com_1 com_2 com_3 com_4 com_1 com_2
## Levels: com_1 com_2 com_3 com_4
col.id.factor.mr[a_mostrar]
## [1] 1 1 1 1 2 2
## 40 Levels: 1 2 5 8 9 10 11 12 14 15 16 17 19 21 22 24 27 28 29 30 31 33 35 36 37 39 40 41 43 ... 59

post_frd_durbin_Jamovi <- durbinAllPairsTest(y = col.v.depend.mr ,
                                             groups = col.v.factor.mr,
                                             blocks = col.id.factor.mr,
                                             p.adjust.method = "none")

post_frd_durbin_holm <- durbinAllPairsTest(y = col.v.depend.mr,
                                            groups = col.v.factor.mr,
                                            blocks = col.id.factor.mr,
                                            p.adjust.method = "holm")
post_frd_durbin_Jamovi
##       com_1   com_2  com_3 
## com_2 0.0325  -      -     
## com_3 2.1e-07 0.0011 -     
## com_4 2.1e-07 0.0011 1.0000

post_frd_durbin_holm
##       com_1   com_2  com_3 
## com_2 0.0651  -      -     
## com_3 1.3e-06 0.0043 -     
## com_4 1.3e-06 0.0043 1.0000

Vemos que, efectivamente, el test que lanzaba Jamovi ofrece unos valores de significación más bajos, resultado de no haber controlado el error FWER, por lo que no son fiables. Además, el test de Durbin está diseñado para el caso de bloques incompletos (incomplete block design) por lo que mejor aplicamos un post hoc como el de Nemenyi o el test más actual de Eisinga (2017), usando la función frdAllPairsExactTest:

post_frd_eisinga <- frdAllPairsExactTest(y = col.v.depend.mr ,
                                         groups = col.v.factor.mr,
                                         blocks = col.id.factor.mr,
                                         p.adjust.method = "holm")

post_frd_eisinga
##       com_1   com_2 com_3
## com_2 0.182   -     -    
## com_3 4.5e-05 0.032 -    
## com_4 4.5e-05 0.032 1.000

Podemos apreciar que hay dos grupos, el formado por {com_1, com_2} y otro formado por {com_3, com_4}. Tanto com_1 como com_2 son marcas con una menor preferencia por parte de los clientes que cualquiera de las marcas com_3 y com_4.

Post hoc genérico

El post hoc usando un test genérico con penalización lo hacemos de una forma análoga al caso paramétrico. En ese caso, habíamos usado la fórmula básica sin término de error y la función pairwise_t_test. Ahora usamos una fórmula análoga y la función pairwise_wilcox_test que aplica el test no paramétrico para 2 muestras pareadas de Wilcoxon (Wilcoxon signed rank test) a todos los pares. Debemos pasar como parámetro paired = TRUE. Por defecto usa la penalización de Holm:

formula.mr.sin.error = as.formula(paste0(v.depend.mr, "~", v.factor.mr))
formula.mr.sin.error
## Valoracion ~ Marca

post_pairwise_wilcox_test_mr = pairwise_wilcox_test(formula.mr.sin.error, 
                                     data = datos.mr.long,
                                     paired = TRUE)  
GraficoCajasPost(datos.mr.long, v.depend.mr, v.factor.mr, 
                 es_param =  FALSE, es_mr = TRUE, post_pairwise_wilcox_test_mr)

Nota. También se puede usar la función wilcox_test. Esta función realiza automáticamente las comparaciones múltiples cuando ve que el factor tiene más de un valor.

De nuevo, podemos apreciar que se distinguen los mismos grupos con medias distintas: el formado por {com_1, com_2} y el formado por {com_3, com_4}.

En cualquier caso, para poder aplicar cada una de las comparaciones múltiples con el Wilcoxon signed rank test deberíamos confirmar que las distribuciones de todas las diferencias son simétricas. Para ello, podemos construir todas las posibles combinaciones de las columnas con CombPairs y realizar un plot de las densidades de dichas diferencias:

comb = CombPairs(columnas.mr)
comb
##      X1    X2
## 1 com_1 com_2
## 2 com_1 com_3
## 3 com_1 com_4
## 4 com_2 com_3
## 5 com_2 com_4
## 6 com_3 com_4

par(mfrow=c(3, 2))
sapply(c(1:nrow(comb)) , 
       function(x) plot (density(datos.mr.wide[comb[x,1]][,1] -  datos.mr.wide[comb[x,2]][,1]) , main = ""))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL

Es aceptable admitir que todas las anteriores distribuciones son simétricas, por lo que el post-hoc genérico también encuentra significativa la diferencia entre los grupos {com_1, com_2} y {com_3, com_4}.

Post hoc con control

También podemos seleccionar un control, fijando un grupo y comparando el resto con él. Para ello, podemos usar la función frdManyOneExactTest que implementa el post hoc con control de Eisinga (también está disponible el de Nemenyi a través de la función frdManyOneNemenyiTest)

IMPORTANTE: El control se establece según el orden lexicográfico entre los nombres de las columnas.

post_frd_eisinga_control = frdManyOneExactTest(y = col.v.depend.mr ,
                                               groups = col.v.factor.mr,
                                               blocks = col.id.factor.mr,
                                               p.adjust.method = "holm")
post_frd_eisinga_control
##       com_1  
## com_2 0.091  
## com_3 2.3e-05
## com_4 2.3e-05

Podemos ver que, al haber fijado un grupo de control y disminuir el número de comparaciones, hemos rebajado el p-value de la comparación entre com_1 y com_2 de 0.182 a 0.091, aunque no lo suficiente como para que sea significativo.

2.3 Test para comparar resultados de varios clasificadores

Vamos a aplicar lo visto anteriormente al conjunto de datos de clasificación. El fichero RData lo puede descargar desde este enlace:

Classification Algorithms (RData)

Vamos a suponer que queremos comparar el algoritmo 3 con el resto. Tal y como hemos visto en la sección anterior, debemos garantizar que el grupo de control tenga el nombre en la posición más baja según un orden lexicográfico. Por tanto, lo que vamos a hacer es anteponerle un símbolo _:

load("ClassificationAlgorithms.RData")
datos.mr.wide = ClassificationAlgorithms
nombres.columnas = colnames(datos.mr.wide)  
v.id.mr = nombres.columnas[1]
v.depend.mr = "TasaAcierto"
v.factor.mr = "Algoritmo"
grupo.control.mr = "Alg3"
indice.control = match(grupo.control.mr, nombres.columnas)
colnames(datos.mr.wide)[indice.control] = paste0("_", grupo.control.mr)
nombres.columnas = colnames(datos.mr.wide)

Para realizar la transformación al formato long le asignamos a columnas.mr todas las columnas menos la primera (la de índice 1) ya que contiene el nombre del dataset y no es una medida repetida. También vamos a suprimir del análisis el algoritmo 1 (columna con índice 2) ya que presenta un comportamiento sistemáticamente peor que el resto y no queremos que la razón de que el test de Friedman resulte significativo sea porque hemos incluido dicho algoritmo en el estudio (que por supuesto va a tener una mediana mucho menor que el resto)

columnas.mr = nombres.columnas[-c(1,2)]

Realizamos la transformación de los datos al formato long y establecemos los nombres de las variables:

datos.mr.long = pivot_longer(datos.mr.wide, 
                             cols = any_of(columnas.mr), 
                             names_to = v.factor.mr,
                             values_to = v.depend.mr)

head(datos.mr.long)
## # A tibble: 6 × 4
##   Dataset  Alg1 Algoritmo TasaAcierto
##   <chr>   <dbl> <chr>           <dbl>
## 1 aud      25.3 Alg2             76  
## 2 aud      25.3 _Alg3            75  
## 3 aud      25.3 Alg4             69.6
## 4 aud      25.3 Alg5             71  
## 5 aud      25.3 Alg6             70.9
## 6 aud      25.3 Alg7             57.7

indice.col.id.mr = which(colnames(datos.mr.long) == v.id.mr)
indice.col.factor.mr = which(colnames(datos.mr.long) == v.factor.mr)
datos.mr.long[,indice.col.id.mr] = as.factor(as.data.frame(datos.mr.long)[,indice.col.id.mr])
datos.mr.long[,indice.col.factor.mr] =as.factor(as.data.frame(datos.mr.long)[,indice.col.factor.mr])

col.v.depend.mr = as.data.frame(datos.mr.long)[,v.depend.mr]
col.v.factor.mr = as.data.frame(datos.mr.long)[,v.factor.mr]
col.id.factor.mr = as.data.frame(datos.mr.long)[,v.id.mr]

Comprobamos la normalidad de los residuos, una vez que hemos eliminado el Algoritmo 1. Recuerde que teníamos que establecer la fórmula formula.anova.mr.error y llamar a la función ResiduosMR:

formula.anova.mr.error = as.formula(paste0(v.depend.mr, "~",
                                    v.factor.mr, "+ Error(" ,
                                    v.id.mr , "/" , v.factor.mr , ")"))
formula.anova.mr.error
## TasaAcierto ~ Algoritmo + Error(Dataset/Algoritmo)

residuos.mr = ResiduosMR(datos.mr.long, v.depend.mr, v.factor.mr,
                         v.id.mr, columnas.mr)
ggqqplot(residuos.mr)

Realmente, los residuos no se alejan demasiado de una distribución normal. Sin embargo, si lanzamos el test de Shapiro, obtenemos que se rechaza la normalidad:

shapiro.test(residuos.mr)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos.mr
## W = 0.93373, p-value = 1.821e-06

Por tanto, si queremos ser rigurosos, en vez de lanzar un ANOVA paramétrico de medidas repetidas, debemos lanzar un test no paramétrico (de Friedman). Por lo tanto, debemos establecer la fórmula pertinente para ejecutar la función friedman_test:

formula.anova.friedman_test = as.formula(paste0(v.depend.mr, "~", 
                                                v.factor.mr , "|" , 
                                                v.id.mr))
formula.anova.friedman_test
## TasaAcierto ~ Algoritmo | Dataset

test_friedman = friedman_test (formula.anova.friedman_test, 
                               data = datos.mr.long)
GraficoCajasTest(datos.mr.long, v.depend.mr, v.factor.mr,
                 es_param = FALSE, es_mr = TRUE ,
                 test_friedman)

Podemos apreciar que el test resulta significativo. Aceptamos por tanto que hay al menos un algoritmo con un comportamiento distinto. Lanzamos el post hoc, aplicando el test de Eisinga, usando la función frdAllPairsExactTest de PMCMRPlus (puede probar también con el test de Nemenyi usando la función frdAllPairsNemenyiTest):

post_frd_eisinga <- frdAllPairsExactTest(y = col.v.depend.mr ,
                                         groups = col.v.factor.mr,
                                         blocks = col.id.factor.mr)
post_frd_eisinga
##      _Alg3 Alg2  Alg4  Alg5  Alg6 
## Alg2 0.032 -     -     -     -    
## Alg4 0.073 1.000 -     -     -    
## Alg5 0.023 1.000 1.000 -     -    
## Alg6 0.023 1.000 1.000 1.000 -    
## Alg7 1.000 0.220 0.391 0.185 0.185

Tenga en cuenta que el algoritmo 3 aparece el primero porque le habíamos cambiado el nombre, pero el test se ha aplicado comparando todos con todos. Podemos observar que ninguna de las comparaciones múltiples resulta significativa, por lo que el post-hoc no ha sido capaz de detectar el algoritmo que es distinto al resto (se ha producido un error de tipo II ya que el test de Friedman sí era significativo). Sin embargo, si sólo estamos interesados en comparar el algoritmo 3 con el resto, podemos fijarlo como control. Como estamos comparando 6 algoritmos, en vez de realizar 15 tests simultáneos (6x5/2), pasaremos a realizar únicamente 5 comparaciones por lo que la penalización del nivel de significación será bastante menor. Para ello, usamos el post hoc de Eisinga con control "_Alg3" (recuerde que el grupo de control es el que aparece primero en un orden alfabético). Nos sale lo siguiente:

post_frd_eisinga_control = frdManyOneExactTest(y = col.v.depend.mr ,
                                               groups = col.v.factor.mr,
                                               blocks = col.id.factor.mr,
                                               p.adjust.method = "holm")
post_frd_eisinga_control
##      _Alg3 
## Alg2 0.0075
## Alg4 0.0121
## Alg5 0.0075
## Alg6 0.0075
## Alg7 0.4702

por lo que concluimos que el algoritmo 3 no se puede diferenciar del 7 pero sí del resto. Finalmente, vamos a ver la distribución de las variables de medidas repetidas (los algoritmos):

ggdensity(as.data.frame(datos.mr.long), v.depend.mr, facet.by = v.factor.mr)

Los algoritmos 2, 3, 4 y 7 tienen distribuciones similares, pero algo distintas a los algoritmos 5 y 6. Por lo tanto, no podemos extraer conclusiones con rigor estadístico sobre las medianas y el test de Friedman sería un test de dominancia estocástica.

Guion de Prácticas de la parte de Introducción a la Estadística del curso Introducción a la Ciencia de Datos

© Juan Carlos Cubero

04 octubre, 2024

1 Tests de Hipótesis con Jamovi

1.1 Tests para muestras independientes

1.1.1 Prueba para dos muestras independientes (test paramétrico)

1.1.2 Test no paramétrico para comparar dos muestras independientes

1.1.3 Anova para k muestras independientes (test paramétrico)

1.1.3.1 Test conjunto (omnibus)

1.1.3.2 Post hoc

1.1.4 Anova para k muestras independientes (test no paramétrico)

1.1.4.1 Test conjunto (omnibus)

1.1.4.2 Post hoc

1.2 Tests para muestras dependientes o de medidas repetidas

1.2.1 Prueba T (t-Student) para dos muestras pareadas (medidas repetidas)

1.2.2 Test no paramétrico para comparar dos muestras dependientes (medidas repetidas)

1.2.3 Anova para k muestras dependientes (test paramétrico de medidas repetidas)

1.2.3.1 Test conjunto (omnibus)

1.2.3.2 Post hoc

1.2.4 Anova para k muestras dependientes (test no paramétrico de medidas repetidas)

1.2.4.1 Test conjunto (omnibus)

1.2.4.2 Post hoc

1.3 Test para comparar resultados de varios clasificadores

2 Tests de Hipótesis con R

2.1 Tests para muestras independientes

2.1.1 Prueba para dos muestras independientes (test paramétrico)

2.1.2 Test no paramétrico para comparar dos muestras independientes

2.1.3 Anova para k muestras independientes (test paramétrico)

2.1.3.1 Test conjunto (omnibus)

2.1.3.2 Post hoc

Post hoc específico

Post hoc genérico

Post hoc con control

2.1.4 Anova para k muestras independientes (test no paramétrico)

2.1.4.1 Test conjunto (omnibus)

2.1.4.2 Post hoc

Post hoc específico

Post hoc genérico

Post hoc con control

2.2 Tests para muestras dependientes o de medidas repetidas

2.2.1 Prueba T (t-Student) para dos muestras pareadas (medidas repetidas)

2.2.2 Test no paramétrico para comparar dos muestras dependientes (medidas repetidas)

2.2.3 Anova para k muestras dependientes (test paramétrico de medidas repetidas)

2.2.3.1 Test conjunto (omnibus)

2.2.3.2 Post hoc

Post hoc específico

Post hoc genérico

Post hoc con control

2.2.4 Anova para k muestras dependientes (test no paramétrico de medidas repetidas)

2.2.4.1 Test conjunto (omnibus)

2.2.4.2 Post hoc

Post hoc específico

Post hoc genérico

Post hoc con control

2.3 Test para comparar resultados de varios clasificadores