Légère incohérence entre la fonction R intégrée de Kruskal-Wallis et le calcul manuel

9

Je suis confus par ce qui suit, et je n'ai pas pu trouver la réponse ailleurs.

J'essaie d'apprendre R tout en faisant quelques statistiques, et, comme exercice, j'essaie de revérifier les résultats des fonctions R intégrées en les faisant également `` à la main '', pour ainsi dire, dans R. Cependant , pour le test de Kruskal-Wallis, je continue d'obtenir des résultats différents, et je ne peux pas comprendre pourquoi.

Par exemple, je regarde les données suivantes distribuées dans un exercice

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)

Et je veux analyser l'activité par groupe. Je lance d'abord un test de Kruskal-Wallis en utilisant la fonction R intégrée

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

Ce qui renvoie .H=8.9056

Pour revérifier, j'essaie de faire la même chose «à la main» dans R, avec le code (sans doute sans défense) suivant

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

Ce qui est censé refléter la formule suivante:

H=12N(N+1)i=1g(Ri2ni)3(N+1)

NgniiRii

H=8.499H

J'ai essayé de chercher à comprendre ce que je fais mal ou ne comprends pas, mais en vain. Quelqu'un peut-il m'aider à comprendre pourquoi la kruskal.testfonction intégrée renvoie une valeur différente de celle que j'obtiens en énonçant les choses?

MSR
la source

Réponses:

12

kruskal.testapplique une correction pour les liens comme décrit dans cet article Wikipedia (point 4):

1i=1G(ti3ti)N3N

Suite de votre code:

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

Vous pouvez découvrir ce que fait la fonction R en étudiant attentivement le code que vous pouvez voir en utilisant getAnywhere(kruskal.test.default).

Roland
la source
4
@MichaelChernick Non, ce n'est pas le cas. Le fait est que OP a appris une simplification du test qui ne devrait être utilisée qu'en l'absence de liens.
Roland
4
@MichaelChernick Je ne dis pas que cela ne rentrerait pas chez Stack Overflow. Mais je dirais que cela convient aussi bien à CV. De toute évidence, il aurait été utile que OP ait non seulement partagé son code, mais aussi les formules qu'il utilise.
Roland
3
@Michael Le statut de ce fil est un appel facile: il est carrément de notre ressort car il cherche à comprendre un test statistique.
whuber
2
Modifié pour inclure la formule reflétée dans le code. J'aurais dû penser à le faire la première fois. Mes excuses.
MSR
3
Voir aussi la fonction de Hmiscpackage R spearman2qui utilise les midranks pour les égalités et un Ftest pour obtenir Kruskal-Wallis. Je pense que c'est plus précis que certaines méthodes.
Frank Harrell