Análisis no paramétrico
library(tidyverse)
library(coin)
library(rstatix)
library(ggpubr)Cuando usar estimación no-paramétrica
La estimación no paramétrica es usada para sustituir estimaciones basadas en la suposición de que las muestras de datos que estamos analizando fueron tomadas de poblaciones que responden a una distribución normal, es decir las muestras presentan una distribución aproximada a la normal o normal, y se usa la media como medida central de muestras poblaciones para estimar sus diferencias. Además otra suposición importante es la homocedasticidad (igualdad de varianzas en las poblaciones).
Se sabe que cuando alguna de estas dos suposiciones, normalidad y homocedasticidad, no se cumplen, esto puede resultar en conclusiones erróneas de nuestras estimaciones. Un posible origen de este problema es resultado de una escala “débil” al medir nuestros datos. En estos casos se usan estimaciones que no están basadas en distribuciones de probabilidad, y las medidas centrales se basan en la mediana. Varias de estas técnicas están basadas en la categorización de los datos y estimación de estadísticas relacionadas con la mediana. Como comentario, en el caso de los algoritmos de R y en lo general se debe tener un mínimo de datos, en este caso por lo menos 6 para que las funciones no marquen un error.
Prueba de signo Wilcoxon categorizada para una muestra.
Esta prueba de Wilcoxon se una para determinar si una muestra no sigue una distribución no paramétrica o fue medida en una escala débil y se quiere determinar si la mediana de la población toma algún valor hipotético o es igual a algún estándar. Para usar esta prueba su supone que los datos están distribuidos cerca de la mediana.
Ejemplo C13_S03_1 Researchers wished to know if instruction in personal care and grooming would improve the appearance of mentally challenged girls. In a school for the mentally challenged, 10 girls selected at random received special instruction in personal care and grooming. Two weeks after completion of the course of instruction, the girls were interviewed by a nurse and a social worker who assigned each girl a score based on her general appearance. The investigators believed that the scores achieved the level of an ordinal scale. They felt that although a score of, say, 8 represented a better appearance than a score of 6, they were unwilling to say that the difference between scores of 6 and 8 was equal to the difference between, say, scores of 8 and 10; or that the difference between scores of 6 and 8 represented twice as much improvement as the difference between scores of 5 and 6. The scores are shown in the following Table. We wish to know if we can conclude that the median score of the population from which we assume this sample to have been drawn is different from 5.
# Tabla de calificaciones de apariencia de 10 niñas en escuela para suetos con capacidades especiales
Apariencia <- tibble(
Sujeto = c(1,2,3,4,5,6,7,8,9,10),
Calif = c(4,5,8,8,9,6,10,7,6,6)
)
Apariencia# A tibble: 10 × 2
Sujeto Calif
<dbl> <dbl>
1 1 4
2 2 5
3 3 8
4 4 8
5 5 9
6 6 6
7 7 10
8 8 7
9 9 6
10 10 6# Estimamos la mediana y el rango intercuartil
Apariencia %>% get_summary_stats(Calif, type = "median_iqr")# A tibble: 1 × 4
variable n median iqr
<chr> <dbl> <dbl> <dbl>
1 Calif 10 6.5 2# Gráfica con los datos
ggboxplot(
Apariencia$Calif, width = 0.5, add = c("mean", "jitter"),
ylab = "Calificación", xlab = FALSE
)
# La prueba de Wilcoxon se corre usando rstatix, que usa la librería coin
calif.test <- Apariencia %>% wilcox_test(Calif ~ 1, mu = 5, alternative = "two.sided")
calif.test# A tibble: 1 × 6
.y. group1 group2 n statistic p
* <chr> <chr> <chr> <int> <dbl> <dbl>
1 Calif 1 null model 10 42.5 0.0197# El tamaño del efecto se obtiene con:
Apariencia %>% wilcox_effsize(Calif ~ 1, mu = 5)# A tibble: 1 × 6
.y. group1 group2 effsize n magnitude
* <chr> <chr> <chr> <dbl> <int> <ord>
1 Calif 1 null model 0.764 10 large Por lo que podemos concluír que la mediana no es igual a \(5\), y que el “efecto” estadístico para llegar a esta conclusión es “grande”, esto último estimado por la función de wilcox_effsize.
Prueba Mann-Whitney
La prueba no paramétrica que sustituye a la prueba t para la comparación de dos poblaciones, tiene varios nombres; Mann-Whitney U, Wilcoxon-Mann-Whitney o prueba Wilcoxon para dos muestas. En R la llamada de la función es a través de la prueba de Wilcoxon wilcox_test.
Para usar esta prueba se deben cumplir las condiciones básicas de muestreo aleatorio, para obtener las muestras poblacionales que se van a usar, y además
1. La escala de medición de las muestras por lo menos debe ser ordinal.
2. La variable de análisis debe ser continua.
3. Si las poblaciones son diferentes solo son diferentes con respecto a sus medianas.
Cuando estas suposiciones se cumplen se puede probar la hipótesis nula \(H_0\) que supone que las poblaciones tienen medianas iguales, en contra de tres hipótesis alternativas \(H_A\) posibles: (1) Las poblaciones tiene medianas distintas (opción “two.sided”), (2) La mediana de la “primera población” es mayor (opción “greater”), (3) La mediana de la “primera población” es menor (opción “less).
Prueba no paramétrica para comparar dos poblaciones (Prueba U Mann-Whitney).
Esta prueba nos permite comparar dos poblaciones por sus medianas, se supone que los datos están cercanos en cada población a sus medianas y la cantidad de datos de cada muestra pueden ser distintas. Si las muestras resultan diferentes solo podemos concluír que son diferentes en sus medianas.
Ejemplo C13.S06.01 A researcher designed an experiment to assess the effects of prolonged inhalation of cadmium oxide. Fifteen laboratory animals served as experimental subjects while 10 similar animals served as controls. The variable of interest was hemoglobin level following the experiment. The results are shown in EXA_C13_S06_01. We wish to know if we can conclude that prolonged inhalation of cadmium oxide reduces hemoglobin level.
Exa13.6.1 <- read_csv(file="~/Dropbox/GitHub/ProbEstad/DataSets/ch13_all/EXA_C13_S06_01.csv", show_col_types = FALSE)
Exa13.6.1 # tibble de X y Y de los niveles de hemoglobina.# A tibble: 15 × 2
X Y
<dbl> <dbl>
1 14.4 17.4
2 14.2 16.2
3 13.8 17.1
4 16.5 17.5
5 14.1 15
6 16.6 16
7 15.9 16.9
8 15.6 15
9 14.1 16.3
10 15.3 16.8
11 15.7 NA
12 16.7 NA
13 13.7 NA
14 15.3 NA
15 14 NA Ca_O <- Exa13.6.1 %>%
rename( Exposed = X, Control = Y)
# Se alargan los datos
Ca_O_long <- Ca_O %>%
pivot_longer(cols = c("Exposed", "Control") , names_to = "Group", values_to = "HemoglobinLevel", values_drop_na = TRUE)
Ca_O_long <- Ca_O_long %>%
mutate( Group = Group %>% fct_relevel("Control", "Exposed"))
# para ver las medianas de cada muestra
Ca_O_long %>%
group_by(Group) %>%
get_summary_stats(HemoglobinLevel, type = "median_iqr")# A tibble: 2 × 5
Group variable n median iqr
<fct> <chr> <dbl> <dbl> <dbl>
1 Control HemoglobinLevel 10 16.6 1
2 Exposed HemoglobinLevel 15 15.3 1.7# Box plot con la nube de puntos
HemoPlot <- ggboxplot(
Ca_O_long, x = "Group", y = "HemoglobinLevel",
ylab = "Hemoglobin level in blood", xlab = "Groups", add = "jitter"
)
# se muestra la gráfica
HemoPlot
Ca_O_stat.test <- Ca_O_long %>%
wilcox_test(HemoglobinLevel ~ Group) %>%
add_significance()
# Imprimimos la tabla para ver el resultado de la prueba
Ca_O_stat.test# A tibble: 1 × 8
.y. group1 group2 n1 n2 statistic p p.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr>
1 HemoglobinLevel Control Exposed 10 15 125 0.00601 ** # Se estima el tamaño del afecto después de tomar o desechar la Ho
Ca_O_long %>% wilcox_effsize(HemoglobinLevel ~ Group)# A tibble: 1 × 7
.y. group1 group2 effsize n1 n2 magnitude
* <chr> <chr> <chr> <dbl> <int> <int> <ord>
1 HemoglobinLevel Control Exposed 0.555 10 15 large # Gráfica con resultado de la prueba de cajas y puntos con tamaño de efecto y tipo de corrección
Ca_O_stat.test <- Ca_O_stat.test %>% add_xy_position(x = "Group")
HemoPlot +
stat_pvalue_manual(Ca_O_stat.test, tip.length = 0) +
labs(subtitle = get_test_label(Ca_O_stat.test, detailed = TRUE))
Prueba de Wilcoxon para dos muestras pareadas.
Cuando los datos a analizar están pareados, es decir dos observaciones por medida (o sujeto) un buen sustituto para una prueba t cuando los datos cumplen con las suposiciones de la prueba t, es la prueba de signo pareada de Wilcoxon.
# Ejemplo C13.S03.2 A dental research team wished to know if teaching people how to brush their teeth would be beneficial. Twelve pairs of patients seen in a dental clinic were obtained by carefully matching on such factors as age, sex, intelligence, and initial oral hygiene scores. One member of each pair received instruction on how to brush his or her teeth and on other oral hygiene matters. Six months later, all 24 subjects were examined and assigned an oral hygiene score by a dental hygienist unaware of which subjects had received the instruction. A low score indicates a high level of oral hygiene. The results are shown in EXA_C13_S03_02.csv. \(X_i\) are instructed pairs and \(Y_i\) are the uninstructed pairs.
Exa13.3.2 <- read_csv(file="~/Dropbox/GitHub/ProbEstad/DataSets/ch13_all/EXA_C13_S03_02.csv", show_col_types = FALSE)
Exa13.3.2# A tibble: 12 × 3
PAIR X Y
<dbl> <dbl> <dbl>
1 1 1.5 2
2 2 2 2
3 3 3.5 4
4 4 3 2.5
5 5 3.5 4
6 6 2.5 3
7 7 2 3.5
8 8 1.5 3
9 9 1.5 2.5
10 10 2 2.5
11 11 3 2.5
12 12 2 2.5Hygiene <- Exa13.3.2 %>%
rename( Instructed = X, Uninstructed = Y)
# Ahora para alargar los datos, solo dejamos fuera PAIR
Hygiene_long <- Hygiene %>%
pivot_longer(cols = c("Instructed", "Uninstructed") , names_to = "Group", values_to = "HygieneScore")
# Hacer factores a los pares de sujetos y los grupos
Hygiene_long <- Hygiene_long %>%
mutate( PAIR = factor(PAIR), Group = factor(Group) )
# Vamos a ordenar de acuerdo a Uninstructed para ver el efecto contra instrucción debe disminuir -> less
Hygiene_long <- Hygiene_long %>%
mutate( Group = Group %>% fct_relevel("Uninstructed", "Instructed") )
# En este caso por la forma en que está pareada la muestra una gráfica de cajas
HyPlot <- ggboxplot(
Hygiene_long, x = "Group", y = "HygieneScore", add = "jitter",
ylab = "Hygiene Score", xlab = "Groups", color = "blue"
)
# mostramos la gráfica
HyPlot
Hygiene_stat.test <- Hygiene_long %>%
wilcox_test(HygieneScore ~ Group, paired = TRUE, alternative = "greater") %>%
add_significance()
Hygiene_stat.test# A tibble: 1 × 8
.y. group1 group2 n1 n2 statistic p p.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr>
1 HygieneScore Uninstructed Instructed 12 12 57 0.0145 * Hygiene_long %>%
wilcox_effsize(HygieneScore ~ Group, paired = TRUE)# A tibble: 1 × 7
.y. group1 group2 effsize n1 n2 magnitude
* <chr> <chr> <chr> <dbl> <int> <int> <ord>
1 HygieneScore Uninstructed Instructed 0.645 12 12 large Hygiene_stat.test <- Hygiene_stat.test %>% add_xy_position(x = "Group")
HyPlot +
stat_pvalue_manual(Hygiene_stat.test, tip.length = 0) +
labs(subtitle = get_test_label(Hygiene_stat.test, detailed= TRUE))
Ejercicio 13.3.2Determining the effects of grapefruit juice on pharmacokinetics of oral digoxin (a drug often prescribed for heart ailments) was the goal of a study by Parker et al. Seven healthy nonsmoking volunteers participated in the study. Subjects took digoxin with water for 2 weeks, no digoxin for 2 weeks, and digoxin with grapefruit juice for 2 weeks. The average peak plasma digoxin concentration (Cmax) when subjects took digoxin with water is given in the first column of the table in EXR_C13_S03_02.csv. The second column gives the Cmax concentration when subjects took digoxin with grapefruit juice. May we conclude on the basis of these data that the Cmax concentration is higher when digoxin is taken with grapefruit juice? Let \(\alpha = .5\).
Exr13.3.2 <- tibble(
Subject = c(1, 2, 3, 4, 5, 6, 7),
H2O = c(2.34, 2.46, 1.87, 3.09, 5.59, 4.05, 6.21),
Juice = c(3.03, 3.46, 1.97, 3.81, 3.07, 2.62, 3.44)
)
Exr13.3.2# A tibble: 7 × 3
Subject H2O Juice
<dbl> <dbl> <dbl>
1 1 2.34 3.03
2 2 2.46 3.46
3 3 1.87 1.97
4 4 3.09 3.81
5 5 5.59 3.07
6 6 4.05 2.62
7 7 6.21 3.44Pharma_long <- Exr13.3.2 %>%
pivot_longer(cols = c("H2O", "Juice") , names_to = "Group", values_to = "dioxinConcen")
Pharma_long <- Pharma_long %>%
mutate( Group = Group %>% fct_relevel("H2O", "Juice") )
# Estimamos la gráfica con las direcciones de la concentración
PharmaPlot <- ggpaired(
Pharma_long, x = "Group", y = "dioxinConcen",
order = c("H2O", "Juice"),
ylab = "Dioxin Concentration in Blood", xlab = "Groups", color = "blue"
)
# mostramos la gráfica
PharmaPlot
Pharma_stat.test <- Pharma_long %>%
wilcox_test(dioxinConcen ~ Group, paired = TRUE, alternative = "less") %>%
add_significance()
Pharma_stat.test# A tibble: 1 × 8
.y. group1 group2 n1 n2 statistic p p.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr>
1 dioxinConcen H2O Juice 7 7 18 0.766 ns Pharma_long %>%
wilcox_effsize(dioxinConcen ~ Group, paired = TRUE)# A tibble: 1 × 7
.y. group1 group2 effsize n1 n2 magnitude
* <chr> <chr> <chr> <dbl> <int> <int> <ord>
1 dioxinConcen H2O Juice 0.256 7 7 small Como podemos ver del valor de \(p\) no podemos desechar la hipótesis nula, y la conclusión es que no son diferentes las poblaciones de donde se extrajeron las muestras. No hay diferencia entre el uso de agua y jugo de toronja.
Prueba Kruskal Wallis
El sustituto no paramétrico de la prueba ANOVA de una vía, es la prueba categórica de Kruskal-Wallis, el proceso se aplica a \(n_1, n_2, \dots, n_k\) observaciones de \(k\) muestras que combinadas son una serie ordenada de mayor a menor, de tamaño \(n\). Posteriormente las observaciones se les categoriza; se les asigna un rango a partir de \(1\) que se le asigna a la observación más pequeña, hasta observación \(n\), a la que se le asigna el rango mayor. Cuando una o más observaciones se repiten se les asigna el promedio del rango a cada una. Posteriormente, los rangos asignados a cada una de las \(k\) observaciones se suman por separado en k-esimas sumas de rangos. Y se calcula la estadística de prueba: \[ H = \frac{12}{n(n+1)} \sum_{j-1}^{k} \frac{R_{j}^2}{n_j} - 3(n+1) \]
con \(k\) número de uestras, \(n_j\) número de observaciones en la j-esima muestra, \(n\) número total de observaciones combinadas, \(R_j\) la suma de rangos de las observaciones de la j-esima muestra.
La hipótesis nula es que las medidas centrales del todas las muestras están alineadas (son “iguales”).
Ejemplo 13.8.1 In a study of pulmonary effects on guinea pigs, Lacroix et al. exposed ovalbumin (OA)-sensitized guinea pigs to regular air, benzaldehyde, or acetaldehyde. At the end of exposure, the guinea pigs were anesthetized and allergic responses were assessed in bronchoalveolar lavage (BAL). One of the outcome variables examined was the count of eosinophil cells, a type of white blood cell that can increase with allergies. The table in the problem gives the eosinophil cell count (× 106) for the three treatment groups.
OAsensit <- tibble(
Air = c(12.22, 28.44, 28.13, 38.69, 54.91),
Benzaldehyde = c(3.86, 4.05, 6.47, 21.12, 3.33),
Acetaldehyde = c(54.36, 27.87, 66.81, 46.27, 30.19) )
OAsensit_long <- OAsensit %>%
pivot_longer(cols = c("Air", "Benzaldehyde", "Acetaldehyde") , names_to = "Group", values_to = "EosinophilCount")
OAsensit_long <- OAsensit_long %>%
mutate( Group = Group %>% fct_relevel("Air", "Benzaldehyde", "Acetaldehyde") )
# Estimación de las medianas y sus cuantiles
OAsensit_long %>%
group_by(Group) %>%
get_summary_stats(EosinophilCount, type = "median_iqr")# A tibble: 3 × 5
Group variable n median iqr
<fct> <chr> <dbl> <dbl> <dbl>
1 Air EosinophilCount 5 28.4 10.6
2 Benzaldehyde EosinophilCount 5 4.05 2.61
3 Acetaldehyde EosinophilCount 5 46.3 24.2 # Una gráfica de barras
ggboxplot(OAsensit_long, x = "Group", y = "EosinophilCount", add = "jitter", fill = "Group")
OAsensit_stat.test <- OAsensit_long %>%
kruskal_test(EosinophilCount ~ Group)
OAsensit_stat.test# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <int> <dbl> <chr>
1 EosinophilCount 15 9.14 2 0.0104 Kruskal-Wallis# y estimación del tamaño del efecto la etha
OAsensit_long %>% kruskal_effsize(EosinophilCount ~ Group)# A tibble: 1 × 5
.y. n effsize method magnitude
* <chr> <int> <dbl> <chr> <ord>
1 EosinophilCount 15 0.595 eta2[H] large Una vez que se determina si hay un efecto, que no todas las medianas son iguales, se debe determinar en donde está el efecto (cual o cuales medianas). Del resultado de la estimación de Kruskal-Wallis se sabe que hay un efecto entre los grupos pero no se sabe entre que grupos hay diferencias. Por lo general una prueba de Kruskal-Wallis significativa se sigue de una prueba de Dunn para estimar diferencias por pares y determinar entre que grupos está la diferencia. Es posible también usar una prueba de Wilcoxon entre pares y hacer una corrección por medidas múltiples.
# Comparación por pares de Dunn
pwc <- OAsensit_long %>%
dunn_test(EosinophilCount ~ Group, p.adjust.method = "bonferroni")
pwc# A tibble: 3 × 9
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
1 Eosinophi… Air Benzald… 5 5 -2.19 0.0284 0.0851 ns
2 Eosinophi… Air Acetald… 5 5 0.707 0.480 1 ns
3 Eosinophi… Benzald… Acetald… 5 5 2.90 0.00374 0.0112 * # o una prueba de Wilcoxon por pares y corrección
pwc2 <- OAsensit_long %>%
wilcox_test(EosinophilCount ~ Group, p.adjust.method = "bonferroni")
pwc2# A tibble: 3 × 9
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
1 EosinophilCount Air Benzal… 5 5 24 0.016 0.048 *
2 EosinophilCount Air Acetal… 5 5 8 0.421 1 ns
3 EosinophilCount Benzal… Acetal… 5 5 0 0.008 0.024 * # Grafica con error y subtítulo con valores de la estimación
pwc <- pwc %>% add_xy_position(x = "Group")
ggboxplot(OAsensit_long, x = "Group", y = "EosinophilCount", fill = "Group") +
stat_pvalue_manual(pwc, hide.ns = TRUE) +
labs(
subtitle = get_test_label(OAsensit_stat.test, detailed = TRUE),
caption = get_pwc_label(pwc)
)
Ejemplo con más de tres muestras y muestras de datos grandes. En el caso de que las muestras tengan varios datos (muestras “grandes”) hay que ajustar la estimación. Este tipo de estimación es ajustado en forma automática en las funciones de R.
Ejemplo 13.8.2 Table in EXA_C13_S08_02.csv shows the net book value of equipment capital per bed for a sample of hospitals from each of five types of hospitals. We wish to determine, by means of the Kruskal–Wallis test, if we can conclude that the average net book value of equipment capital per bed differs among the five types of hospitals. The ranks of the 41 values, along with the sum of ranks for each sample, are shown in the file table.
Exa13.8.2 <- read_csv(file="~/Dropbox/GitHub/ProbEstad/DataSets/ch13_all/EXA_C13_S08_02.csv", show_col_types = FALSE)
Exa13.8.2# A tibble: 10 × 5
A B C D E
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1735 5260 2790 3475 6090
2 1520 4455 2400 3115 6000
3 1476 4480 2655 3050 5894
4 1688 4325 2500 3125 5705
5 1702 5075 2755 3275 6050
6 2667 5225 2592 3300 6150
7 1575 4613 2601 2730 5110
8 1602 4887 1648 NA NA
9 1530 NA 1700 NA NA
10 1698 NA NA NA NAHosp_long <- Exa13.8.2 %>%
pivot_longer(cols = A:E, names_to = "Type", values_to = "Value", values_drop_na = TRUE)
Hosp_long <- Hosp_long %>%
mutate(Type = Type %>% fct_relevel("A","B","C","D","E"))
Hosp_long %>% group_by(Type) %>%
get_summary_stats(Value, type = "median_iqr")# A tibble: 5 × 5
Type variable n median iqr
<fct> <chr> <dbl> <dbl> <dbl>
1 A Value 10 1645 160.
2 B Value 8 4750 639.
3 C Value 9 2592 255
4 D Value 7 3125 205
5 E Value 7 6000 270.ggboxplot(Hosp_long, x = "Type", y = "Value", add = "jitter", fill = "Type" )
Hosp_stat.test <- Hosp_long %>%
kruskal_test(Value ~ Type)
Hosp_stat.test# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <int> <dbl> <chr>
1 Value 41 36.4 4 0.00000024 Kruskal-Wallis# Comparación Dunn por pares para determinar en donde estan las diferencias
Hosp_pwc <- Hosp_long %>%
dunn_test(Value ~ Type, p.adjust.method = "bonferroni")
Hosp_pwc# A tibble: 10 × 9
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
1 Value A B 10 8 4.21 0.0000250 0.000250 ***
2 Value A C 10 9 1.27 0.205 1 ns
3 Value A D 10 7 2.70 0.00702 0.0702 ns
4 Value A E 10 7 5.24 0.000000163 0.00000163 ****
5 Value B C 8 9 -2.92 0.00355 0.0355 *
6 Value B D 8 7 -1.30 0.195 1 ns
7 Value B E 8 7 1.12 0.261 1 ns
8 Value C D 9 7 1.48 0.139 1 ns
9 Value C E 9 7 3.97 0.0000734 0.000734 ***
10 Value D E 7 7 2.34 0.0191 0.191 ns Hosp_pwc <- Hosp_pwc %>% add_xy_position(x = "Type")
ggboxplot(Hosp_long, x = "Type", y = "Value", add = "jitter", fill = "Type") +
stat_pvalue_manual(Hosp_pwc, hide.ns = TRUE) +
labs(
subtitle = get_test_label(Hosp_stat.test, detailed = TRUE),
caption = get_pwc_label(Hosp_pwc)
)
Prueba de Friedman
El sustituto no paramétrico de la prueba ANOVA de dos vías es la prueba de Friedman, en la que se calcula la estadística \(\chi_r^2\) para estimar la probabilidad. En este caso se tienen \(n\) bloques y \(k\) tratamientos. Para estimar la estadística se suman las categorías de cada variable estas son las \(R_j\) la la estadística de Friedman \(\chi_r^2\) se estima con
\[ \chi_r^2 = \frac{12}{nk(k+1)} \sum_{j=1}^k (R_j)^2 - 3n(k+1) \] Para el caso de la prueba de Friedman, la W de Kendall se usa para estimar el tamaño del efecto. El coeficiente W de Kendall toma valores entre \(0\) (indicando no relación alguna) y \(1\) (indicando una relación perfecta), que sigue la indicación de Cohen de valores entre 0.1 - 0.3 como efecto pequeño, 0.3 - 0.5 para efecto moderado y mayor a 0.5 para un efecto grande, en R se usa friedman_effectsize. Para hacer el post hoc con comparaciones múltiples corregidas por Bonferroni se usa la prueba de signo categorizada de Wilcoxon para identificar que tratamiento es diferente.
Example 13.9.1 A physical therapist conducted a study to compare three models of low-volt electrical stimulators. Nine other physical therapists were asked to rank the stimulators in order of preference. A rank of 1 indicates first preference. The results are shown in EXA_C13_S09_01.csv. We wish to know if we can conclude that the models are not preferred equally.
Exa13.9.1 <- read_csv(file="~/Dropbox/GitHub/ProbEstad/DataSets/ch13_all/EXA_C13_S09_01.csv", show_col_types = FALSE)
Exa13.9.1# A tibble: 9 × 4
THERAPIST A B C
<dbl> <dbl> <dbl> <dbl>
1 1 2 3 1
2 2 2 3 1
3 3 2 3 1
4 4 1 3 2
5 5 3 2 1
6 6 1 2 3
7 7 2 3 1
8 8 1 3 2
9 9 1 3 2ElectStim_long <- Exa13.9.1 %>%
pivot_longer(cols = A:C, names_to = "Stimulator", values_to = "Preference")
ElectStim_long <- ElectStim_long %>%
mutate(Stimulator = Stimulator %>% fct_relevel("A","B","C"))
ElectStim_long %>%
group_by(Stimulator) %>%
get_summary_stats(Preference, type = "median_iqr")# A tibble: 3 × 5
Stimulator variable n median iqr
<fct> <chr> <dbl> <dbl> <dbl>
1 A Preference 9 2 1
2 B Preference 9 3 0
3 C Preference 9 1 1ElectStim_stat.test <- ElectStim_long %>%
friedman_test(Preference ~ Stimulator | THERAPIST)
ElectStim_stat.test# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <dbl> <dbl> <chr>
1 Preference 9 8.22 2 0.0164 Friedman testElectStim_long %>% friedman_effsize(Preference ~ Stimulator | THERAPIST)# A tibble: 1 × 5
.y. n effsize method magnitude
* <chr> <int> <dbl> <chr> <ord>
1 Preference 9 0.457 Kendall W moderate pwc <- ElectStim_long %>%
wilcox_test(Preference ~ Stimulator, paired = TRUE, p.adjust.method = "bonferroni")
pwc# A tibble: 3 × 9
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
1 Preference A B 9 9 3.5 0.023 0.07 ns
2 Preference A C 9 9 24.5 0.851 1 ns
3 Preference B C 9 9 42 0.021 0.062 ns Ejercicio 13.9.1 The table in EXR_C13_S09_01 shows the scores made by nine randomly selected student nurses on final examinations in three subject areas. Test the null hypothesis that student nurses constituting the population from which the sample was drawn perform equally well in all three subject areas against the alternative hypothesis that they perform better in, at least, one area. Let \(\alpha = .05\).
Exr13.9.1 <- read_csv(file="~/Dropbox/GitHub/ProbEstad/DataSets/ch13_all/EXR_C13_S09_01.csv", show_col_types = FALSE)
Exr13.9.1# A tibble: 9 × 4
SUBJ FUND PHYS ANAT
<dbl> <dbl> <dbl> <dbl>
1 1 98 95 77
2 2 95 71 79
3 3 76 80 91
4 4 95 81 84
5 5 83 77 80
6 6 99 70 93
7 7 82 80 87
8 8 75 72 81
9 9 88 81 83Nurses_long <- Exr13.9.1 %>%
pivot_longer(cols = c("FUND", "PHYS", "ANAT"), names_to = "Areas", values_to = "Grades")
Nurses_long <- Nurses_long %>% mutate(Areas = Areas %>% fct_relevel("FUND", "PHYS", "ANAT"))
Nurses_long %>%
group_by(Areas) %>%
get_summary_stats(Grades, type = "median_iqr")# A tibble: 3 × 5
Areas variable n median iqr
<fct> <chr> <dbl> <dbl> <dbl>
1 FUND Grades 9 88 13
2 PHYS Grades 9 80 9
3 ANAT Grades 9 83 7Nurses_stat.test <- Nurses_long %>%
friedman_test(Grades ~ Areas | SUBJ)
Nurses_stat.test# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <dbl> <dbl> <chr>
1 Grades 9 8.67 2 0.0131 Friedman testNurses_long %>% friedman_effsize(Grades ~ Areas | SUBJ)# A tibble: 1 × 5
.y. n effsize method magnitude
* <chr> <int> <dbl> <chr> <ord>
1 Grades 9 0.481 Kendall W moderate Nurses_stat.test_EffSize <- Nurses_long %>%
wilcox_test(Grades ~ Areas, paired = TRUE, p.adjust.method = "bonferroni")
Nurses_stat.test_EffSize# A tibble: 3 × 9
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
1 Grades FUND PHYS 9 9 41 0.033 0.098 ns
2 Grades FUND ANAT 9 9 31 0.342 1 ns
3 Grades PHYS ANAT 9 9 8 0.097 0.291 ns Se puede concluir que de la prueba Friedman, hay diferencias en las calificaciones entre las materias de estudio de las enfermeras, sin embargo al estimar la prueba para el tamaño del efecto, esté es moderado (0.481) ya que no llega al .5. Además al estimar la prueba de Wilcoxon por pares y corrección de Bonferroni las diferencias entre las áreas no resultan significativas. La tendencia está entre la diferencia de Fundamentales y Fisiología.