Nieznaczna niespójność między wbudowaną funkcją R Kruskala-Wallisa a obliczeniami ręcznymi

Jestem zdezorientowany następującymi informacjami i nie byłem w stanie znaleźć odpowiedzi gdzie indziej.

Próbuję nauczyć się języka R podczas wykonywania statystyk i, jako ćwiczenie, próbuję dwukrotnie sprawdzić wyniki wbudowanych funkcji R, wykonując je również „ręcznie”, jak to było w R. Jednak , dla testu Kruskala-Wallisa ciągle otrzymuję różne wyniki i nie mogę zrozumieć, dlaczego.

Na przykład patrzę na następujące dane przekazane w ćwiczeniu

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)

Chcę analizować aktywność według grup. Najpierw uruchamiam test Kruskala-Wallisa przy użyciu wbudowanej funkcji R.

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

Który zwraca $H = 8.9056$.

Aby sprawdzić dwukrotnie, próbuję zrobić to samo „ręcznie” w R, używając następującego (bez wątpienia bezradnego) kodu

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

Które ma odzwierciedlać następującą formułę:

Gdzie $N$ oznacza całkowitą liczbę obserwacji, $g$ to liczba grup, $n_i$ to liczba obserwacji w $i$grupa oraz $R_i$ jest sumą stopni $i$grupa.

A teraz rozumiem $H = 8.499$, co, dodatkowo do mojego zamieszania, jest również odpowiedzią na dane ćwiczenie. Próbowałem tego dla kilku różnych zestawów danych i mam tendencję do uzyskiwania nieco wyższej wartości$H$ za pomocą wbudowanej funkcji.

Próbowałem dowiedzieć się, co robię źle lub nie rozumiem, ale bezskutecznie. Czy ktoś może mi pomóc zrozumieć, dlaczego wbudowana kruskal.testfunkcja zwraca inną wartość niż ta, którą otrzymuję, wypowiadając rzeczy?

kruskal.teststosuje korektę powiązań, jak opisano w tym artykule w Wikipedii (punkt 4):

Korekcję powiązań, jeśli stosuje się wzór skrótu opisany w poprzednim punkcie, można wykonać, dzieląc H przez $1 - \frac{\sum_{i=1}^G (t_i^3 - t_i)}{N^3-N}$, ...

Kontynuując od kodu:

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

Możesz dowiedzieć się, co robi funkcja R, uważnie studiując kod, którego możesz użyć getAnywhere(kruskal.test.default).