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łę:

H = \frac{12}{N (N + 1)} \sum_{i = 1}^{g} (\frac{R_{i}^{2}}{n_{i}}) - 3 (N + 1)

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

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?

r kruskal-wallis ties

— MSR
źródło

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).

— Roland
źródło

@MichaelChernick Nie, to nie jest. Chodzi o to, że OP został nauczony uproszczenia testu, który powinien być stosowany tylko wtedy, gdy nie ma żadnych powiązań.

— Roland

@MichaelChernick Nie twierdzę, że nie zmieściłby się w przepełnieniu stosu. Ale twierdziłbym, że równie dobrze pasuje do CV. Oczywiście byłoby pomocne, gdyby OP dzielił nie tylko swój kod, ale także formuły, których używają.

— Roland

@Michael Status tego wątku jest łatwy do sprawdzenia: jest on całkowicie zgodny z naszym zakresem, ponieważ ma on na celu zrozumienie testu statystycznego.

— whuber

Edytowane w celu uwzględnienia formuły odzwierciedlonej w kodzie. Powinienem był to zrobić za pierwszym razem. Przeprosiny.

— MSR

Zobacz także funkcję Hmiscpakietu R , spearman2która wykorzystuje midranki do wiązania i Ftest na uzyskanie Kruskala-Wallisa. Myślę, że jest to dokładniejsze niż niektóre metody.

— Frank Harrell,