Interpretación de la prueba de inmersión de Hartigans

Question

Me gustaría encontrar una manera de cuantificar la intensidad de la bimodalidad de algunas distribuciones que obtuve empíricamente. Por lo que he leído, todavía hay cierto debate sobre la forma de cuantificar la bimodalidad. He optado por utilizar el test de inmersión de Hartigans que parece ser el único disponible en R (paper original : http://www.stat.washington.edu/wxs/Stat593-s03/Literature/hartigan85a.pdf ). La prueba de inmersión de Hartigans se define como : "La prueba de inmersión mide la multimodalidad en una muestra por la diferencia máxima, sobre todos los puntos de la muestra, entre la función de distribución empírica, y la función de distribución unimodal que minimiza esa diferencia máxima" .

Me gustaría entender completamente cómo debo interpretar estas estadísticas antes de utilizarlas. Esperaba que la prueba de inmersión aumentara si la distribución es multimodal (ya que se define como "la máxima diferencia con respecto a la distribución unimodal"). Pero : se puede leer en la página de la wikipedia sobre la distribución multimodal que "Los valores inferiores a 0,05 indican una bimodalidad significativa y los valores superiores a 0,05 pero inferiores a 0,10 sugieren una bimodalidad con significación marginal". . Dicha afirmación proviene de este papel (Fig. 2). Según este documento, el índice de la prueba de inmersión se acerca a 0 cuando la distribución es bimodal. Esto me confunde.

Para interpretar correctamente la prueba de inmersión de Hartigans construí algunas distribuciones (el código original es de aquí ) y aumenté el valor de exp(mu2) (llamado "Intensidad de la bimodularidad" a partir de ahora - Edición : Debería haberlo llamado "Intensidad de la bimodalidad". ) para obtener la bimodalidad. En el primer gráfico, puedes ver algunos ejemplos de distribuciones. Luego estimé el índice diptest (segundo gráfico) y el valor p (tercer gráfico) asociado (paquete diptest ) a esas diferentes distribuciones simuladas. El código R utilizado está al final de mi post.

Lo que muestro aquí es que el índice de la prueba de inmersión es alto y el valor P es bajo cuando las distribuciones son bimodales. Lo cual es contrario a lo que se puede leer en Internet.

No soy un experto en estadística, por lo que apenas he entendido el artículo de Hartigans. Me gustaría recibir algunos comentarios sobre la forma correcta de interpretar la prueba de inmersión de Hartigans. ¿Estoy equivocado en alguna parte?

Gracias a todos. Saludos,

T.A.

Ejemplo de distribución simulada :

Índice de prueba de inmersión de Hartigan asociado :

Prueba de inmersión de Hartigan valor p.asociado :

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

Answer 1

El Sr. Freeman (autor del papel que te comenté) me dijo que en realidad sólo se fijaba en el valor P de la prueba de inmersión. Esta confusión proviene de su frase :
"Los valores del HDS van de 0 a 1, con valores inferiores a 0,05 que indican una bimodalidad significativa, y valores superiores a 0,05 pero inferiores a 0,10 que sugieren una bimodalidad con significación marginal" . Los valores HDS corresponden al valor P, y no a la estadística de la prueba de inmersión. No estaba claro en el documento.

Mi análisis es bueno: la estadística de la prueba de inmersión aumenta cuando la distribución se desvía de una distribución unimodal.

La prueba de bimodalidad y la prueba de Silverman también se pueden calcular fácilmente en R y hacen bien el trabajo.