Test hipotez dla różnicy w medianach pomiędzy więcej niż dwiema próbkami

12

Pytanie

Wyniki testu trzech grup ludzi są zapisywane jako osobne wektory w R.

set.seed(1)
group1 <- rnorm(100, mean = 75, sd = 10)
group2 <- rnorm(100, mean = 85, sd = 10)
group3 <- rnorm(100, mean = 95, sd = 10)

Chcę wiedzieć, czy istnieje znacząca różnica w medianach między tymi grupami. Wiem, że mogłem przetestować grupę 1 w porównaniu z grupą 2 za pomocą testu Wilcoxona.

wilcox.test(group1, group2)

Jednak porównuje to tylko dwie grupy na raz i chciałbym porównać wszystkie trzy jednocześnie. Chciałbym przeprowadzić test statystyczny, który daje wartość ap na poziomie istotności 0,05. Czy ktoś mógłby pomóc?

Edycja # 1 - Mediana testu Mooda

Po sugerowanej odpowiedzi użytkownika Hibernating spróbowałem testu mediany Mooda.

median.test <- function(x, y){
    z <- c(x, y)
    g <- rep(1:2, c(length(x), length(y)))
    m <- median(z)
    fisher.test(z < m, g)$p.value
}

median.test(group1, group2)

Jednak takie podejście pozwala mi przetestować istotną różnicę między medianami tylko dwóch grup jednocześnie. Nie jestem pewien, jak go użyć, aby porównać mediany wszystkich trzech jednocześnie.

Edycja # 2 - Test Kruskala-Wallisa

Sugerowana odpowiedź użytkownika dmartin wydaje się być mniej więcej tym, czego potrzebuję i pozwala mi przetestować wszystkie trzy grupy jednocześnie.

kruskal.test(list(group1, group2, group3))

Edytuj # 3

Użytkownik Greg Snow z pomocą zauważa w swojej odpowiedzi, że test Kruskala-Wallisa jest odpowiedni, o ile zawiera on surowe założenia, które czynią go również testem środków.

Alexander
źródło
Na tej stronie pojawiło się już wiele podobnych pytań. Proszę poszukać median test. Moja własna odpowiedź / komentarze jest tutaj .
ttnphns
Jeśli chodzi o porównanie median wszystkich trzech jednocześnie, zobacz moją edycję nieco zmodyfikowanego kodu R.
Hibernacja

Odpowiedzi:

4

Można również zastosować test Kruskala-Wallisa, ponieważ jest to nieparametryczna ANOVA. Ponadto często uważa się go za silniejszy niż test mediany Mooda . Można go zaimplementować w języku R za pomocą funkcji kruskal.test w pakiecie statystyk w języku R.

Aby odpowiedzieć na Twoją edycję, interpretacja KW jest podobna do jednostronnej ANOVA. Znacząca wartość p odpowiada odrzuconemu zeru, że wszystkie trzy średnie są równe. Musisz użyć testu kontrolnego (ponownie, podobnie jak ANOVA), aby odpowiedzieć na pytania dotyczące określonych grup. Zwykle wynika to z określonych pytań badawczych. Wystarczy spojrzeć na parametry symulacji, wszystkie trzy grupy powinny się znacznie różnić od siebie, jeśli wykonasz test kontrolny (ponieważ wszystkie są 1 SD z wyjątkiem N = 100).

dmartin
źródło
1
Aby wyjaśnić kilka rzeczy. 1) Kruskal-Wallis nie jest testem median, chyba że rozkłady obserwacji w grupach spełniają pewne założenia. Jeśli naprawdę chcesz porównać mediany, może to nie być odpowiedni test. Najlepiej wybrać test, który faktycznie testuje hipotezę, którą chcesz przetestować. 2) Kruskal-Wallis nie jest „ANOVA”. Oznacza to, że nie jest to analiza wariancji. 3) W tej odpowiedzi wzmianka o „środkach” jest nieprawidłowa.
Sal Mangiafico
10

Po pierwsze, test Wilcoxona (lub test Manna-Whitneya) nie jest testem median (chyba że przyjmujesz bardzo surowe założenia, które również czynią z niego test środków). Dla porównania więcej niż 2 grup test Wilcoxona może prowadzić do paradoksalnych wyników (patrz Kości Efrona ).

Ponieważ test Wilcoxona jest tylko specjalnym przypadkiem testu permutacji i jesteś szczególnie zainteresowany medianami, proponuję test permutacji na medianach.

Najpierw wybierz miarę różnicy, coś w rodzaju największej z 3 median minus najmniejsza z 3 (lub wariancji 3 median lub MAD itp.).

Teraz oblicz swoje statystyki dla oryginalnych danych.

połącz wszystkie dane w jeden zestaw, a następnie losowo podziel wartości na 3 grupy

te same rozmiary co oryginał i oblicz te same statystyki.

powtarzaj wiele razy (jak 9998)

Porównaj porównanie statystyki z rzeczywistych danych z rozkładem wszystkich statystyk dla twojego testu.

Greg Snow
źródło
Powiedzmy, że jestem gotów przyjąć ścisłe założenia niezbędne do testu Wilcoxa, które uczynią go również testem środków. Czy wymagałoby to zmiany kodu R, który napisałem powyżej? Czy można to zrobić również w przypadku testu Kruskala-Wallisa?
Alexander
1
@Alexander, jeśli chcesz przyjąć te założenia, kod R jest w porządku, a Kruskal Wallis również będzie w porządku. Ale jeśli zechcesz przyjąć te założenia, to t.testi aovprawdopodobnie również będzie dobrze.
Greg Snow
+1. Jeśli mówisz o, Wilcoxon sum-rank testczy nie miałbyś nic przeciwko konwersji „Wilcox” na tę nazwę?
ttnphns
@GregSnow +1 za zdobyte punkty ... ale zakładam, że przez „Wilcox” masz na myśli test nazwany imieniem Franka Wilcoxona. (To zamieszanie jest niestety spotęgowane przez R, który - wprowadzając w błąd - wywołuje odpowiedni test wilcox.test). Mógłbyś edytować?
Glen_b
8

Mediana testu Mooda to test nieparametryczny, który służy do testowania równości median z dwóch lub więcej populacji. Zobacz tutaj część R swojego pytania. Zobacz także powiązane pytanie tutaj . Również stąd :

Test mediany Mooda jest najłatwiejszy do wykonania ręcznie: oblicz ogólną medianę (wszystkich danych) i policz, ile wartości jest powyżej i poniżej mediany w każdej grupie. Jeśli grupy są prawie takie same, obserwacje powinny wynosić około 50-50 powyżej i poniżej ogólnej mediany w każdej grupie ... Liczby poniżej mediany i powyżej mediany ... tworzą tabelę dwukierunkową, która jest następnie analizowany za pomocą testu chi-kwadrat. Mediana testu Mooda jest bardzo podobna do testu znaku uogólnionego na dwie lub więcej grup.

Edycja: w przypadku trzech grup możesz rozważyć tę prostą generalizację kodu R, do którego linkowałem:

median.test2 <- function(x, y, z) {
  a <- c(x, y, z)
  g <- rep(1:3, c(length(x), length(y), length(z)))
  m <- median(a)
  fisher.test(a < m, g)$p.value
}
Hibernacja
źródło
1
+1 za nazwanie testu. Nie wiedziałem, że test mediany nazywa się również testem Mooda.
ttnphns
+1 Dzięki za pomoc w tym, naprawdę to doceniam!
Alexander
Wiem o kilku implementacjach R. mood.medtestw pakiecie RVAideMemoire wydaje się być zwykłym testem, z tym wyjątkiem, że domyślnie używa on dokładnego testu Fishera dla mniejszych próbek. median_testFunkcji w pakiecie monety może stanowić asymptotycznej testu, lub Monte Carlo.
Sal Mangiafico
0

Wiem, że jest już późno, ale nie mogłem też znaleźć dobrego pakietu dla testu mediany Mooda, więc wziąłem na siebie zadanie stworzenia funkcji w R, która wydaje się załatwić sprawę.

#Mood's median test for a data frame with one column containing data (d),
#and another containing a factor/grouping variable (f)

moods.median = function(d,f) {

    #make a new matrix data frame
    m = cbind(f,d)
    colnames(m) = c("group", "value")


    #get the names of the factors/groups
    facs = unique(f)

    #count the number of factors/groups
    factorN = length(unique(f))


    #Make a 2 by K table that will be saved to the global environment by using "<<-":
    #2 rows (number of values > overall median & number of values <= overall median)
    #K-many columns for each level of the factor
    MoodsMedianTable <<- matrix(NA, nrow = 2, ncol = factorN)

    rownames(MoodsMedianTable) <<- c("> overall median", "<= overall median")
    colnames(MoodsMedianTable) <<- c(facs[1:factorN])
    colnames(MoodsMedianTable) <<- paste("Factor: ",colnames(MoodsMedianTable))


    #get the overall median
    overallmedian = median(d)



    #put the following into the 2 by K table:
    for(j in 1:factorN){ #for each factor level

        g = facs[j] #assign a temporary "group name"


        #count the number of observations in the factor that are greater than
        #the overall median and save it to the table
        MoodsMedianTable[1,j] <<- sum(m[,2][ which(m[,1]==g)] > overallmedian)


        #count the number of observations in the factor that are less than
        # or equal to the overall median and save it to the table
        MoodsMedianTable[2,j] <<- sum(m[,2][ which(m[,1]==g)] <= overallmedian)

    }


    #percent of cells with expected values less than 5
    percLT5 = ((sum(chisq.test(MoodsMedianTable)$expected < 5)) /
        (length(chisq.test(MoodsMedianTable)$expected)))


    #if >20% of cells have expected values less than 5
    #then give chi-squared stat, df, and Fisher's exact p.value
    if (percLT5 > 0.2) {
        return(list(
            "Chi-squared" = chisq.test(MoodsMedianTable)$statistic,
            "df" = chisq.test(MoodsMedianTable)$parameter,
            "Fisher's exact p.value" = fisher.test(MoodsMedianTable)$p.value))

    }


    #if <= 20% of cells have expected values less than 5
    #then give chi-squared stat, df, and chi-squared p.value
    if (percLT5 <= 0.2) {
        return(list(
            "Chi-squared" = chisq.test(MoodsMedianTable)$statistic,
            "df" = chisq.test(MoodsMedianTable)$parameter,
            "Chi-squared p.value" = chisq.test(MoodsMedianTable)$p.value))

    }

}

W przypadku pytania PO najpierw uruchom to, aby utworzyć nową ramkę danych, w której będą przechowywane wartości z trzech wektorów grupy za pomocą dopasowanej zmiennej „group”.

require(reshape2)
df = cbind(group1, group2, group3)
df = melt(df)
colnames(df) = c("observation", "group", "value")

i uruchom funkcję testu mediany Mooda za pomocą moods.median(df$value, df$group)

JRF1111
źródło
Wydaje się, że odpowiedzią był test Kruskala-Wallisa. PO potrzebował rozwiązania z 3 grupami. Wygląda na to, że ttnphns już dostarczył kod R do testu nastroju.
Michael R. Chernick
1
Kod podany przez ttnphns zapewnia tylko wartość ap, ten, który napisałem, podaje również chi kwadrat stat i df i działa dla dowolnej liczby grup. Przeważnie właśnie tu pisałem, ponieważ ten post jest pierwszym, który pojawia się, gdy szukam sposobu wykonania mediany testu mediany
Mooda