Interpretacja testu Hartigansa ' dip test
On 18 listopada, 2020 by adminChciałbym znaleźć sposób na ilościowe określenie intensywności bimodalności niektórych dystrybucje, które otrzymałem empirycznie. Z tego, co przeczytałem, wciąż toczy się dyskusja na temat sposobu kwantyfikacji bimodalności. Zdecydowałem się użyć testu zanurzeniowego Hartigansa, który wydaje się być jedynym dostępnym na R (oryginalny artykuł: http://www.stat.washington.edu/wxs/Stat593-s03/Literature/hartigan85a.pdf ). „Test zanurzeniowy Hartigansa jest definiowany jako: „ Test zanurzeniowy mierzy multimodalność w próbce na podstawie maksymalnej różnicy, we wszystkich punktach próbkowania, między funkcją dystrybucji empirycznej a rozkładem unimodalnym, która minimalizuje tę maksymalną różnicę ”.
Chciałbym dokładnie zrozumieć, jak powinienem interpretować te statystyki przed ich użyciem. Spodziewałem się, że test zapadu wzrośnie, jeśli rozkład jest multimodalny (ponieważ jest definiowany jako „maksymalna różnica w stosunku do rozkładu unimodalnego”). Ale : na stronie wikipedii o dystrybucji multimodalnej możesz przeczytać, że „Wartości mniejsze niż 0,05 wskazują na znaczną bimodalność i wartości większa niż 0,05, ale mniejsza niż 0,10 sugeruje bimodalność z marginalnym znaczeniem. ”. Takie stwierdzenie pochodzi z tego artykułu (ryc. 2). Zgodnie z tym artykułem, wskaźnik zapadu jest bliski 0, gdy rozkład jest bimodalny. To mnie dezorientuje.
Aby poprawnie zinterpretować test dip-Hartigansa, skonstruowałem kilka dystrybucji (oryginalny kod pochodzi z tutaj ) i zwiększyłem wartość exp (mu2) (nazywaną od teraz „Intensywnością bimodularności” – Edycja: powinienem był nazwać to „Intensywnością bimodalności” ), aby uzyskać bimodalność. Na pierwszym wykresie można zobaczyć jakiś przykład dystrybucji. Następnie oszacowałem indeks diptest (drugi wykres) i wartość p (trzeci wykres) związane (pakiet diptest ) z tymi różnymi symulowanymi dystrybucjami. Użyty kod R znajduje się na końcu mój post.
Pokazuję tutaj, że indeks testu zanurzenia jest wysoki, a wartość P jest niska, gdy dystrybucje są bimodalne. Co jest sprzeczne z tym, co możesz przeczytać w internecie.
Nie jestem ekspertem w statystykach, więc ledwo rozumiałem artykuł Hartigansa. Chciałbym otrzymać kilka komentarzy na temat właściwej interpretacji testu dip-test Hartigansa. Czy gdzieś się mylę?
Dziękuję wszystkim. Pozdrawiam,
TA
Przykład symulowanej dystrybucji:
Powiązany indeks testu upadu Hartigana:
Powiązana wartość testu zanurzenia Hartigana:
library(diptest) library(ggplot2) # CONSTANT PARAMETERS sig1 <- log(3) sig2 <- log(3) cpct <- 0.5 N=1000 #CREATING BIMOD DISTRIBUTION bimodalDistFunc <- function (n,cpct, mu1, mu2, sig1, sig2) { y0 <- rlnorm(n,mean=mu1, sd = sig1) y1 <- rlnorm(n,mean=mu2, sd = sig2) flag <- rbinom(n,size=1,prob=cpct) y <- y0*(1 - flag) + y1*flag } #DIP TEST DIP_TEST <- function(bimodalData) { TEST <- dip.test(bimodalData) return(TEST$statistic[[1]]) # return(TEST$p.value[[1]]) to get the p value } DIP_TEST(bimodalData) # SIMULATION exp_mu1 = 1 max_exp_mu2 = 100 intervStep = 100 repPerInt = 10 # single distibutions expMu2Value <- c() bimodalData <- c() mu1 <- log(exp_mu1) mu2 <- log(exp_mu1) bimodalData <- c(bimodalData,log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2))) expMu2Value <- c(expMu2Value,rep(exp_mu1,length(log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2))))) mu1 <- log(exp_mu1) mu2 <- log(max_exp_mu2) bimodalData <- c(bimodalData,log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2))) expMu2Value <- c(expMu2Value,rep(max_exp_mu2,length(log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2))))) mu1 <- log(exp_mu1) mu2 <- log(trunc((max_exp_mu2-exp_mu1)/2+1)) bimodalData <- c(bimodalData,log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2))) expMu2Value <- c(expMu2Value,rep(trunc((max_exp_mu2-exp_mu1)/2+1),length(log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2))))) tableExamples <- data.frame(expMu2Value,bimodalData) tableExamples$expMu2Value <- as.factor(tableExamples$expMu2Value) ExamplePlot <- ggplot(tableExamples)+ geom_histogram(aes(bimodalData),color="white")+ ylab("Count")+ xlab("")+ facet_wrap(~expMu2Value)+ ggtitle("Intensity of bimodularity") # calculation of the dip test index exp_mu2Int = seq(from=exp_mu1,to=max_exp_mu2,length.out=intervStep) expmu2Vec = c() dipStat = c() testDone = c() for(exp_mu2 in exp_mu2Int){ mu1 <- log(exp_mu1) mu2 <- log(exp_mu2) for(rep in 1:repPerInt){ bimodalData <- log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2)) diptestone = DIP_TEST(bimodalData) expmu2Vec = c(expmu2Vec,exp_mu2) dipStat = c(dipStat,diptestone) testDone = c(testDone,"diptest") } } table = data.frame(expmu2Vec,dipStat,testDone) IndexPlot <- ggplot(table)+ geom_point(aes(expmu2Vec,dipStat,color=testDone))+ ylab("Index")+ xlab("Intensity of Bimodularity")+ scale_color_discrete(name="Test") ExamplePlot IndexPlot
Komentarze
- Bardzo szczegółowe pytania dotyczące tematu, który jest tajemniczy według standardów ' statystyki . Oczywiste pierwsze pytanie, zanim jeszcze przejdziemy do interpretacji , brzmi: " Dlaczego potrzebujesz tego testu? Jakie informacje ma on przekazać? " Może dostarczyć dodatkowego kontekstu dla motywacji, które doprowadziły cię do znacznie, znacznie późniejszego problemu interpretacja wyników " testu dip? " Innymi słowy, innymi słowy ' s wygodę w programowaniu R, jaka ścieżka logiczna doprowadziła Cię do " testu dip " w pierwszej kolejności?
- Dziękuję za odpowiedź, Mike. Pracuję ' nad modelem teoretycznym w biologii ewolucyjnej i przeprowadzam analizę wrażliwości. W szczególności zauważam, że zmiana niektórych parametrów modyfikuje rozkład zmiennej wyjściowej z unimodalnej na bimodalną (co jest w rzeczywistości bardzo interesujące). To ' jest powodem, dla którego ' szukam prostej statystyki opisującej wielomodularność rozkładu. Pozwoliłoby mi to skoncentrować analizę wrażliwości na multimodularności.
- Dowiedziałem się, że test dip można łatwo obliczyć w R i że może on ilościowo określić odchylenie od rozkładu unimodalnego. Oczywiście byłbym naprawdę zainteresowany wszelkimi innymi statystykami opisującymi wielomodularność rozkładu.
- Hmmm … dopasowanie kilku skromnych wielomianów może dać " biedny człowiek ' s " podejście do radzenia sobie z krzywizną, którą ' ponownie obserwujesz i może być łatwiejsze do wdrożenia i zinterpretowania niż test Hartigana '. Nie ' nie mówisz, czy twoje problemy obejmują zajmowanie się funkcjami wzrostu.Na przykład w rozwoju człowieka istnieje kilka dobrze znanych " wybić " na trajektorii wzrostu w różnych punktach cyklu życia. Stwierdzono, że modele nieparametryczne lepiej pasują i przybliżają te nieliniowości niż modele parametryczne.
- W kwestiach statystycznych: Jak już powiedziano, test upadu przyjmuje jako odniesienie unimodalność. Nie ' nie sądzę, aby odstępstwa od tego można było interpretować w kategoriach liczby trybów tylko na podstawie wartości P. Uważam, że ' znacznie bardziej przydatne jest interpretowanie liczby modów z kombinacją estymacji gęstości i interpretacji merytorycznej.
Odpowiedź
Panie Freeman (autor artykułu , o którym wam mówiłem) powiedział mi, że tak naprawdę patrzył tylko na wartość P testu zanurzeniowego. To zamieszanie bierze się z jego zdania:
„Wartości HDS mieszczą się w zakresie od 0 do 1 z wartościami mniejszymi niż 0,05 wskazującymi na znaczną bimodalność, a wartościami większymi niż 0,05, ale mniejszymi niż 0,10 sugerującymi bimodalność z marginalnym znaczeniem” . Wartości HDS odpowiadają wartości P, a nie statystykom testu zapadu. W artykule było to niejasne.
Moja analiza jest dobra: statystyki testu zapadu rosną, gdy rozkład jest odbiegający od rozkładu unimodalnego.
Test bimodalności i test Silvermana można również łatwo obliczyć w języku R i dobrze wykonać swoje zadanie.
Komentarze
- Zarejestruj & i połącz swoje konta. Informacje o tym, jak to zrobić, można znaleźć w sekcji Moje konto w naszym centrum pomocy .
Dodaj komentarz