Je voudrais trouver un moyen de quantifier l'intensité de la bimodalité de certaines distributions que j'ai obtenues empiriquement. D'après ce que j'ai lu, il y a encore un débat sur la façon de quantifier la bimodalité. J'ai choisi d'utiliser le test d'immersion de Hartigans qui semble être le seul disponible sur R (document original: http://www.stat.washington.edu/wxs/Stat593-s03/Literature/hartigan85a.pdf ). Le test d'immersion de Hartigans est défini comme suit: "Le test d'immersion mesure la multimodalité dans un échantillon par la différence maximale, sur tous les points d'échantillonnage, entre la fonction de distribution empirique et la fonction de distribution unimodale qui minimise cette différence maximale" .
Je voudrais bien comprendre comment interpréter ces statistiques avant de les utiliser. Je m'attendais à ce que le test d'immersion augmente si la distribution est multimodale (car elle est définie comme "la différence maximale par rapport à la distribution unimodale"). Mais : vous pouvez lire dans la page wikipedia sur la distribution multimodale que "les valeurs inférieures à 0,05 indiquent une bimodalité significative et les valeurs supérieures à 0,05 mais inférieures à 0,10 suggèrent une bimodalité avec une signification marginale." . Une telle déclaration vient de cet article (Fig. 2). Selon cet article, l'indice de test d'immersion est proche de 0 lorsque la distribution est bimodale. Ça me perturbe.
Pour interpréter correctement le dip dip de Hartigans, j'ai construit quelques distributions (le code original est d' ici ) et j'ai augmenté la valeur de exp (mu2) (appelée `` Intensité de bimodularité '' à partir de maintenant - Edit: j'aurais dû l'appeler 'Intensité de bimodalité ' ) pour obtenir la bimodalité. Dans le premier graphique, vous pouvez voir quelques exemples de distributions. Ensuite, j'ai estimé l'indice de diptest (deuxième graphique) et la valeur p (troisième graphe) associée (package diptest ) à ces différentes distributions simulées. Le code R utilisé est à la fin de mon message.
Ce que je montre ici, c'est que l'indice de test d'immersion est élevé et que la valeur P est faible lorsque les distibutions sont bimodales. Ce qui est contraire à ce que vous pouvez lire sur Internet.
Je ne suis pas un expert en statistique, de sorte que je comprenais à peine l'article de Hartigans. J'aimerais obtenir quelques commentaires sur la bonne façon d'interpréter le test d'immersion de Hartigans. Ai-je tort quelque part?
Merci à tous. Cordialement,
TA
Exemple de distribution simulée:
Indice de test d'immersion de Hartigan associé:
Test d'immersion de Hartigan valeur p associée:
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
Réponses:
M. Freeman (auteur de l' article dont je vous ai parlé) m'a dit qu'il ne regardait en fait que la valeur du test d'immersion. Cette confusion vient de sa phrase:
"Les valeurs HDS vont de 0 à 1 avec des valeurs inférieures à 0,05 indiquant une bimodalité significative, et des valeurs supérieures à 0,05 mais inférieures à 0,10 suggérant une bimodalité avec une signification marginale" . Les valeurs HDS correspondent à la valeur P et non aux statistiques du test d'immersion. Ce n'était pas clair dans le document.
Mon analyse est bonne: les statistiques du dip test augmentent lorsque la distribution s'écarte d'une distribution unimodale.
Le test de bimodalité et le test de Silverman peuvent également être calculés facilement dans R et font bien le travail.
la source