Ligera inconsistencia entre la función R incorporada de Kruskal-Wallis y el cálculo manual

Question

Estoy confundido por lo siguiente, y no he podido encontrar la respuesta en otra parte.

Estoy tratando de aprender R mientras hago algunas estadísticas, y, como ejercicio, trato de comprobar los resultados de las funciones R incorporadas haciendo también estas "a mano", por así decirlo, en R. Sin embargo, para la prueba de Kruskal-Wallis sigo obteniendo resultados diferentes, y no puedo averiguar por qué.

Por ejemplo, estoy mirando los siguientes datos repartidos en un ejercicio

activity <- c(2, 4, 3, 2, 3, 3, 4, 0, 4, 3, 4, 0, 0, 1, 3, 1, 2, 0, 3, 1, 0, 3, 4, 0, 1, 2, 2, 2, 3, 2) 
group <- c(rep("A", 11), rep("B", 10), rep("C", 9))
group <- factor(group)
data.raw <- data.frame(activity, group)

Y quiero analizar la actividad por grupo. Primero hago una prueba de Kruskal-Wallis usando la función R incorporada

kruskal.test(activity ~ group, data = data.raw)

Que regresa $H = 8.9056$ .

Para comprobarlo, intento hacer lo mismo "a mano" en la R, con el siguiente código (sin duda indefenso)

rank <- rank(activity)
data.rank <- data.frame(rank, group)
rank.sum <- aggregate(rank ~ group, data = data.rank, sum)

x <- rank.sum[1,2]^2 / 11 + rank.sum[2,2]^2 / 10 + rank.sum[3,2]^2 / 9
H <- (12 / (length(activity) * (length(activity) + 1))) * x - 3 * (length(activity) + 1)
H

Que se supone que refleja la siguiente fórmula:

$$ H = \frac {12}{N(N+1)} \sum_ {i = 1}^g \left ( \frac {R^2_i}{n_i} \right ) - 3(N + 1)$$

Donde $N$ es el número total de observaciones, $g$ es el número de grupos, $n_i$ es el número de observaciones en el $i$ el grupo, y $R_i$ es la suma de los rangos de la $i$ el grupo.

Y ahora tengo $H = 8.499$ que, para mi confusión, es también la respuesta dada para el ejercicio en cuestión. He intentado esto para un par de conjuntos de datos diferentes, y tiendo a obtener un valor ligeramente más alto para $H$ usando la función incorporada.

He tratado de buscar para averiguar lo que hago mal o no entiendo, pero no ha servido de nada. ¿Alguien puede ayudarme a entender por qué el edificio kruskal.test devuelve un valor diferente del que obtengo al deletrear las cosas?

Answer 1

kruskal.test aplica una corrección para los empates como se describe en este artículo de Wikipedia (punto 4):

La corrección de los empates si se utiliza la fórmula abreviada descrita en el punto anterior puede hacerse dividiendo H por $1 - \frac{\sum_{i=1}^G (t_i^3 - t_i)}{N^3-N}$ , ...

Continuando con su código:

TIES <- table(activity)
H / (1 - sum(TIES^3 - TIES)/(length(activity)^3 - length(activity)))
#[1] 8.9056

Puede averiguar lo que hace la función de R estudiando detenidamente el código, que puede ver utilizando getAnywhere(kruskal.test.default) .