Interpretace Hartigans ' dip test
On 18 listopadu, 2020 by adminChtěl bych najít způsob, jak kvantifikovat intenzitu bimodality některých distribuce jsem získal empiricky. Z toho, co jsem četl, stále existuje nějaká debata o způsobu kvantifikace bimodality. Rozhodl jsem se použít Hartiganův „dip test, který se zdá být jediným dostupným na R (původní papír: http://www.stat.washington.edu/wxs/Stat593-s03/Literature/hartigan85a.pdf ). Hartigansův „dip test je definován jako: “ Dip test měří multimodalitu ve vzorku maximálním rozdílem ve všech bodech vzorkování mezi empirickou distribuční funkcí a unimodální distribuční funkcí, která tento maximální rozdíl minimalizuje „.
Chtěl bych úplně pochopit, jak mám tuto statistiku interpretovat, než ji použiji. Očekával jsem, že poklesový test se zvýší, pokud je distribuce multimodální (protože je definována jako „maximální rozdíl od unimodální distribuce“). Ale : na stránce wikipedia si můžete přečíst o multimodální distribuci, že „Hodnoty menší než 0,05 označují významnou bimodalitu a hodnoty více než 0,05, ale méně než 0,10 naznačuje bimodalitu s okrajovým významem. “. Takové prohlášení pochází z tohoto článku (obr. 2). Podle tohoto článku je index dip testu téměř 0, když je distribuce bimodální. Matou mě to.
Abych správně interpretoval Hartiganův „dip test, sestrojil jsem několik distribucí (původní kód je z zde ) a zvýšil jsem hodnota exp (mu2) (od nynějška nazývaná „Intenzita bimodularity“ – Upravit: Měl jsem ji nazvat „Intenzita bimodality“ ), abyste získali bimodalitu. V prvním grafu vidíte nějaký příklad distribucí. Poté jsem odhadl index diptestu (druhý graf) a hodnotu p (třetí graf) související (balíček diptest ) s těmito různými simulovanými distribucemi. Použitý kód R je na konci můj příspěvek.
Zde ukážu, že index dip testu je vysoký a Pvalue je nízká, když jsou distribuce bimodální. Což je v rozporu s tím, co si můžete přečíst na internetu.
Nejsem žádný odborník na statistiku, takže jsem sotva rozuměl Hartigansovu práci. Chtěl bych dostat pár komentářů o správném způsobu, jakým bychom měli interpretovat Hartiganův „dip test. Mýlím se někde?
Děkuji všem. S pozdravem,
TA
Příklad simulace distribuce:
Hartiganův index dip testu spojený:
Hartiganův dip test spojený s hodnotou:
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
Komentáře
- Vypracování velmi důkladné otázky týkající se tématu, které je tajemné podle jakýchkoli statistických ' s standardů . Zřejmé první otázky, než se člověk vůbec dostane do interpretace , jsou: " Proč potřebujete tento test? Jaké informace jsou určeny ke komunikaci? " Může poskytnout další kontext pro motivace, které vás vedly k mnohem, mnohem dále navazující otázce interpretace výsledků " dip testu? " Jinými slovy, kromě toho ' s programováním R, jaká cesta logiky vás vedla k " dip testu " na prvním místě?
- Děkuji za odpověď, Miku. ' Pracuji na teoretickém modelu v evoluční biologii a provádím analýzu citlivosti. Zejména pozoruji, že měnící se některé parametry mění distribuci výstupní proměnné z unimodální na bimodální (což je ve skutečnosti velmi zajímavé). Proto ' proto hledám ' jednoduchou statistiku k popisu multimodularity distribuce. Umožnilo by mi to zaměřit analýzu citlivosti na multimodularitu.
- Zjistil jsem, že dip test lze snadno vypočítat v R a že dokáže kvantifikovat odchylku od unimodálního rozdělení. Samozřejmě by mě opravdu zajímala jakákoli další statistika popisující multimodularitu distribuce.
- Hmmm … přizpůsobení několika skromných polynomů by mohlo představovat " chudák ' s " přístup k řešení křivočarosti, kterou ' pozorujete a může být snadněji nasazeno a interpretováno než Hartiganův ' s test. ' neřeknete, zda mezi vaše problémy patří řešení jakýchkoli růstových funkcí.Například v lidském vývoji existuje několik dobře známých " boulí " v trajektorii růstu v různých bodech životního cyklu. Bylo zjištěno, že neparametrické modely lépe pasují a aproximují tyto nelinearity než parametrické modely.
- Ke statistickým problémům: Jak již bylo řečeno, dip test bere jako referenci unimodalitu. ' si nemyslím, že odchylky od toho lze interpretovat z hlediska počtu režimů pouze z P-hodnoty. ' Zjistil jsem, že je mnohem užitečnější interpretovat řadu režimů pomocí kombinace odhadu hustoty a věcné interpretace.
Odpověď
Mr. Freeman (autor článku , o kterém jsem vám řekl) mi řekl, že se ve skutečnosti dívá pouze na Pvalue dip testu. Tento zmatek pochází z jeho věty:
„Hodnoty HDS se pohybují od 0 do 1 s hodnotami menšími než 0,05, což naznačuje významnou bimodalitu, a hodnoty většími než 0,05, ale menšími než 0,10, což naznačuje bimodalitu s okrajovým významem“ . Hodnoty HDS odpovídají Pvalue, nikoli statistikám dip testu. V novinách to bylo nejasné.
Moje analýza je dobrá: statistika dip testu se zvyšuje, když se distribuce odchyluje od unimodální distribuce.
Test bimodality a test Silvermanu lze také snadno vypočítat v prostředí R a dělat práci dobře.
Komentáře
- Zaregistrujte se & a spojte své účty. Informace o tom, jak to provést, najdete v sekci Můj účet v centrum nápovědy .
Napsat komentář