Como testar se minha distribuição é multimodal?

21

Quando plogo um histograma dos meus dados, ele tem dois picos:

histograma

Isso significa uma potencial distribuição multimodal? Eu executei o dip.testem R ( library(diptest)) e a saída é:

D = 0.0275, p-value = 0.7913

Posso concluir que meus dados têm uma distribuição multimodal?

DADOS

10346 13698 13894 19854 28066 26620 27066 16658  9221 13578 11483 10390 11126 13487 
15851 16116 24102 30892 25081 14067 10433 15591  8639 10345 10639 15796 14507 21289 
25444 26149 23612 19671 12447 13535 10667 11255  8442 11546 15958 21058 28088 23827 
30707 19653 12791 13463 11465 12326 12277 12769 18341 19140 24590 28277 22694 15489 
11070 11002 11579  9834  9364 15128 15147 18499 25134 32116 24475 21952 10272 15404 
13079 10633 10761 13714 16073 23335 29822 26800 31489 19780 12238 15318  9646 11786 
10906 13056 17599 22524 25057 28809 27880 19912 12319 18240 11934 10290 11304 16092 
15911 24671 31081 27716 25388 22665 10603 14409 10736  9651 12533 17546 16863 23598 
25867 31774 24216 20448 12548 15129 11687 11581
user1260391
fonte
3
Use mais caixas no seu histograma. Eu sugiro cerca de duas vezes mais
Glen_b -Reinstar Monica
1
Há menção a nove testes diferentes nesta resposta , alguns dos quais podem ser relevantes para sua situação.
Glen_b -Reinstala Monica
1
É provável que este documento seja útil para você, se você ainda não o viu (também neste acompanhamento ) #
303 Eoin

Respostas:

15

O @NickCox apresentou uma estratégia interessante (+1). No entanto, posso considerá-lo de natureza mais exploratória, devido à preocupação que @whuber aponta .

Deixe-me sugerir outra estratégia: você pode ajustar um modelo de mistura finita gaussiana. Observe que isso pressupõe fortemente que seus dados são extraídos de uma ou mais normais verdadeiras. Como @whuber e @NickCox apontam nos comentários, sem uma interpretação substantiva desses dados - apoiada por uma teoria bem estabelecida - para apoiar essa suposição, essa estratégia também deve ser considerada exploratória.

Primeiro, vamos seguir a sugestão de @ Glen_b e examinar seus dados usando o dobro de posições:

insira a descrição da imagem aqui

Ainda vemos dois modos; se alguma coisa, eles aparecem mais claramente aqui. (Observe também que a linha de densidade do kernel deve ser idêntica, mas parece mais espalhada devido ao maior número de posições.)

Agora vamos ajustar um modelo de mistura finita gaussiana. Em R, você pode usar o Mclustpacote para fazer isso:

library(mclust)
x.gmm = Mclust(x)
summary(x.gmm)
# ----------------------------------------------------
# Gaussian finite mixture model fitted by EM algorithm 
# ----------------------------------------------------
#   
# Mclust V (univariate, unequal variance) model with 2 components:
#   
#   log.likelihood   n df       BIC       ICL
#        -1200.874 120  5 -2425.686 -2442.719
# 
# Clustering table:
#  1  2 
# 68 52 

Dois componentes normais otimizam o BIC. Para comparação, podemos forçar um ajuste de um componente e executar um teste de razão de verossimilhança:

x.gmm.1 = Mclust(x, G=1)
logLik(x.gmm.1)
# 'log Lik.' -1226.241 (df=2)
logLik(x.gmm)-logLik(x.gmm.1)
# 'log Lik.' 25.36657 (df=5)
1-pchisq(25.36657, df=3)  # [1] 1.294187e-05

Isso sugere que é extremamente improvável que você encontre dados tão unodais quanto os seus se eles vierem de uma única distribuição normal verdadeira.

Algumas pessoas não se sentem confortáveis ​​usando um teste paramétrico aqui (embora, se as premissas se mantiverem, não conheço nenhum problema). Uma técnica amplamente aplicável é usar o Método de ajuste cruzado de inicialização paramétrica (descrevo o algoritmo aqui ). Podemos tentar aplicá-lo a estes dados:

x.gmm$parameters
# $mean
# 12346.98 23322.06 
# $variance$sigmasq
# [1]  4514863 24582180
x.gmm.1$parameters
# $mean
# [1] 17520.91
# $variance$sigmasq
# [1] 43989870

set.seed(7809)
B = 10000;    x2.d = vector(length=B);    x1.d = vector(length=B)
for(i in 1:B){
  x2      = c(rnorm(68, mean=12346.98, sd=sqrt( 4514863)), 
              rnorm(52, mean=23322.06, sd=sqrt(24582180)) )
  x1      = rnorm( 120, mean=17520.91, sd=sqrt(43989870))
  x2.d[i] = Mclust(x2, G=2)$loglik - Mclust(x2, G=1)$loglik
  x1.d[i] = Mclust(x1, G=2)$loglik - Mclust(x1, G=1)$loglik
}
x2.d = sort(x2.d);  x1.d = sort(x1.d)
summary(x1.d)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# -0.29070 -0.02124  0.41460  0.88760  1.36700 14.01000 
summary(x2.d)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  9.006  23.770  27.500  27.760  31.350  53.500 

insira a descrição da imagem aqui

As estatísticas de resumo e os gráficos de densidade do kernel para as distribuições de amostragem mostram vários recursos interessantes. A probabilidade de log para o modelo de componente único raramente é maior que a do ajuste de dois componentes, mesmo quando o verdadeiro processo de geração de dados tem apenas um componente e, quando é maior, a quantidade é trivial. A idéia de comparar modelos que diferem em sua capacidade de ajustar dados é uma das motivações por trás do PBCM. As duas distribuições de amostragem quase não se sobrepõem; apenas 0,35% x2.dsão inferiores ao valor máximox1.dvalor. Se você selecionasse um modelo de dois componentes se a diferença na probabilidade de log fosse> 9,7, você selecionaria incorretamente o modelo de um componente 0,01% e o modelo de dois componentes 0,02% do tempo. Estes são altamente discrimináveis. Se, por outro lado, você optar por usar o modelo de um componente como hipótese nula, o resultado observado será suficientemente pequeno para não aparecer na distribuição empírica da amostra em 10.000 iterações. Podemos usar a regra 3 (veja aqui ) para colocar um limite superior no valor-p, ou seja, estimamos que seu valor-p seja menor que 0,0003. Ou seja, isso é altamente significativo.

p<.000001p<.001) e os componentes subjacentes, caso existam, também não são garantidos como perfeitamente normais. Se você acha razoável que seus dados possam vir de uma distribuição positiva inclinada, e não normal, esse nível de bimodalidade pode estar dentro do intervalo típico de variação, que é o que suspeito que o teste de mergulho esteja dizendo.

- Reinstate Monica
fonte
1
O problema com essa abordagem é que a alternativa com a qual você está comparando uma mistura gaussiana não é muito razoável. Um mais razoável seria que a distribuição seja algum tipo de distribuição correta, como um Gamma. É quase certo que uma mistura se encaixe em praticamente qualquer conjunto de dados distorcido "significativamente" melhor do que um único gaussiano se ajustará a ele.
whuber
Você está certo, @whuber. Eu tentei fazer esse ponto explicitamente. Não sei como fazer um FMM gama, mas seria melhor.
gung - Restabelece Monica
1
Como isso é exploratório, um pensamento seria tentar transformar a distribuição original em simetria (talvez com uma transformação de Box-Cox deslocada, estimada com robustez a partir de alguns quantis dos dados) e tentar sua abordagem novamente. É claro que você não falaria sobre "significado" em si, mas a análise da probabilidade ainda poderia ser reveladora.
whuber
@ Whuber, fiz isso, mas só mencionei de passagem. (A transformação ideal de Box-Cox é a raiz quadrada inversa.) Você obtém o mesmo resultado, mas os valores de p são (ainda altamente, mas) menos significativos.
gung - Restabelece Monica
3
Eu gosto muito da ideia de que você deveria modelar o que você acha que é o processo de geração. Meu problema é que, mesmo quando as misturas gaussianas funcionam bem, sinto que deve haver uma interpretação substantiva. Se o OP nos dissesse mais sobre quais são os dados, algumas suposições melhores poderiam ser possíveis.
Nick Cox
10

No seguimento das ideias em @ resposta e os comentários de Nick, você pode ver o quão grande as necessidades de largura de banda para ser a apenas alisar o modo secundário:

insira a descrição da imagem aqui

Pegue essa estimativa de densidade do kernel como nulo proximal - a distribuição mais próxima dos dados, mas ainda consistente com a hipótese nula de que é uma amostra de uma população unimodal - e simule a partir dela. Nas amostras simuladas, o modo secundário nem sempre parece tão distinto, e você não precisa ampliar tanto a largura de banda para nivelá-la.

<code> insira a descrição da imagem aqui </code>

A formalização dessa abordagem leva ao teste realizado em Silverman (1981), "Usando estimativas de densidade de núcleo para investigar a modalidade", JRSS B , 43 , 1. O silvermantestpacote de Schwaiger & Holzmann implementa esse teste e também o procedimento de calibração descrito por Hall & York ( 2001), "Na calibração do teste de Silverman para multimodalidade", Statistica Sinica , 11 , p 515, que ajusta o conservadorismo assintótico. A realização do teste em seus dados com uma hipótese nula de unimodalidade resulta em valores de p de 0,08 sem calibração e 0,02 com calibração. Não estou familiarizado o suficiente com o teste de mergulho para adivinhar por que ele pode ser diferente.

Código R:

  # kernel density estimate for x using Sheather-Jones method to estimate b/w:
density(x, kernel="gaussian", bw="SJ") -> dens.SJ
  # tweak b/w until mode just disappears:
density(x, kernel="gaussian", bw=3160) -> prox.null
  # fill matrix with simulated samples from the proximal null:
x.sim <- matrix(NA, nrow=length(x), ncol=10)
for (i in 1:10){
  x.sim[ ,i] <- rnorm(length(x), sample(x, size=length(x), replace=T), prox.null$bw)
}
  # perform Silverman test without Hall-York calibration:
require(silvermantest)
silverman.test(x, k=1, M=10000, adjust=F)
  # perform Silverman test with Hall-York calibration:
silverman.test(x, k=1, M=10000, adjust=T)
Scortchi - Restabelecer Monica
fonte
+1. Interessante! Qual kernel está sendo usado aqui? Pelo que me lembro vagamente, há razões sutis pelas quais os núcleos gaussianos devem ser usados ​​para variantes formais dessa abordagem.
Nick Cox
@ Nick: Kernel gaussiano, mas não me lembro se há uma razão convincente para isso. Cada amostra simulada é redimensionada e há uma correção para um viés conservador do teste original - elaborado por alguém chamado Storey, eu acho.
Scortchi - Restabelecer Monica
@NickCox: Desculpe, não é Storey.
Scortchi - Restabelece Monica
@ Scortchi, alterei um pouco o seu texto e código. Espero que você não se importe. +1. Além disso, você usa o temido operador de atribuição de seta para a direita ?! Oh, a humanidade ...
gung - Reinstate Monica
2
Não é melhor ou pior, realmente, mas a convenção na programação é declarar suas variáveis ​​à esquerda e ter o que lhes é atribuído à direita. Muitas pessoas estão assustadas com-> ; Estou apenas confuso.
gung - Restabelece Monica
7

As coisas para se preocupar incluem:

  1. O tamanho do conjunto de dados. Não é pequeno, nem grande.

  2. A dependência do que você vê na origem do histograma e na largura da lixeira. Com apenas uma opção evidente, você (e nós) não tem idéia de sensibilidade.

  3. A dependência do que você vê no tipo e largura do kernel e quaisquer outras opções feitas na estimativa de densidade. Com apenas uma opção evidente, você (e nós) não tem idéia de sensibilidade.

Em outros lugares, sugeri provisoriamente que a credibilidade dos modos é suportada (mas não estabelecida) por uma interpretação substantiva e pela capacidade de discernir a mesma modalidade em outros conjuntos de dados do mesmo tamanho. (Maior é melhor também ....)

Não podemos comentar sobre nenhum deles aqui. Uma pequena possibilidade de repetibilidade é comparar o que você obtém com amostras de bootstrap do mesmo tamanho. Aqui estão os resultados de um experimento de token usando o Stata, mas o que você vê é arbitrariamente limitado aos padrões do Stata, que são documentados como arrancados do ar . Eu tenho estimativas de densidade para os dados originais e para 24 amostras de inicialização do mesmo.

A indicação (nem mais, nem menos) é o que acho que analistas experientes adivinhariam de qualquer maneira do seu gráfico. O modo esquerdo é altamente repetível e o direito é claramente mais frágil.

Observe que há uma inevitabilidade nisso: como há menos dados próximos do modo direito, eles nem sempre reaparecem em uma amostra de autoinicialização. Mas este também é o ponto chave.

insira a descrição da imagem aqui

Observe que o ponto 3. acima permanece intocado. Mas os resultados estão em algum lugar entre unimodal e bimodal.

Para os interessados, este é o código:

clear 
set scheme s1color 
set seed 2803 

mat data = (10346, 13698, 13894, 19854, 28066, 26620, 27066, 16658, 9221, 13578, 11483, 10390, 11126, 13487, 15851, 16116, 24102, 30892, 25081, 14067, 10433, 15591, 8639, 10345, 10639, 15796, 14507, 21289, 25444, 26149, 23612, 19671, 12447, 13535, 10667, 11255, 8442, 11546, 15958, 21058, 28088, 23827, 30707, 19653, 12791, 13463, 11465, 12326, 12277, 12769, 18341, 19140, 24590, 28277, 22694, 15489, 11070, 11002, 11579, 9834, 9364, 15128, 15147, 18499, 25134, 32116, 24475, 21952, 10272, 15404, 13079, 10633, 10761, 13714, 16073, 23335, 29822, 26800, 31489, 19780, 12238, 15318, 9646, 11786, 10906, 13056, 17599, 22524, 25057, 28809, 27880, 19912, 12319, 18240, 11934, 10290, 11304, 16092, 15911, 24671, 31081, 27716, 25388, 22665, 10603, 14409, 10736, 9651, 12533, 17546, 16863, 23598, 25867, 31774, 24216, 20448, 12548, 15129, 11687, 11581)
set obs `=colsof(data)' 
gen data = data[1,_n] 

gen index = . 

quietly forval j = 1/24 { 
    replace index = ceil(120 * runiform()) 
    gen data`j' = data[index]
    kdensity data`j' , nograph at(data) gen(xx`j' d`j') 
} 

kdensity data, nograph at(data) gen(xx d) 

local xstuff xtitle(data/1000) xla(10000 "10" 20000 "20" 30000 "30") sort 
local ystuff ysc(r(0 .0001)) yla(none) `ystuff'   

local i = 1 
local colour "orange" 
foreach v of var d d? d?? { 
    line `v' data, lc(`colour') `xstuff'  `ystuff' name(g`i', replace) 
    local colour "gs8" 
    local G `G' g`i' 
    local ++i 
} 

graph combine `G' 
Nick Cox
fonte
+1 Gosto da sua abordagem de inicialização: a matriz de gráficos ajuda a todos a entender melhor os dados. Gostaria de saber se esses gráficos podem ser sensíveis à forma como o Stata estima a largura de banda. Suspeito que isso possa resultar em um teste de baixa potência, porque sua estimativa provavelmente se baseia em uma suposição unimodal, levando a uma largura de banda relativamente ampla. Mesmo uma estimativa de largura de banda um pouco mais estreita pode tornar o segundo modo mais proeminente em todas as amostras de inicialização.
whuber
2
@whuber Obrigado! Como sempre, você se concentra infalivelmente nas fraquezas com as quais precisamos nos preocupar, e eu concordo. À medida que a largura de banda do kernel aumenta, a aparência de unimodalidade tende a inevitabilidade. Por outro lado, pequenas larguras de banda geralmente apenas indicam modos espúrios irrepetíveis e / ou triviais. A troca é realmente delicada. Acho que o principal mérito dessa abordagem é a retórica de "Isso é replicável se balançarmos?" Costumo me preocupar com a disposição dos usuários de software de copiar resultados padrão sem reflexão.
Nick Cox
2
Existem abordagens sistemáticas para esse problema com base na modificação progressiva da largura de banda e no rastreamento da aparência e desaparecimento dos modos conforme a largura de banda varia. Em essência, um modo credível persiste e um modo menos credível não. É uma abordagem atraente, mas às vezes seria acionar um construtor de túnel quando uma pá funcionaria. Por exemplo, se você modificar as opções do histograma e o modo secundário desaparecer (ou se mover) com muita facilidade, não acredite.
Nick Cox
2

Identificação do modo não paramétrico LP

Identificação não paramétrica do modo LP (nome do algoritmo LPMode , a ref do artigo é fornecida abaixo)

MaxEntModos [triângulos de cor vermelha na plotagem]: 12783.36 e 24654.28.

L2Modos [triângulos de cor verde na plotagem]: 13054.70 e 24111.61.

Interessante notar as formas modais, especialmente a segunda que mostra considerável distorção (o modelo da Mistura Gaussiana Tradicional provavelmente falhará aqui).

Mukhopadhyay, S. (2016) Identificação do modo em larga escala e ciências orientadas a dados. https://arxiv.org/abs/1509.06428

Deep Mukherjee
fonte
1
Você pode elaborar e fornecer algum contexto para introduzir e explicar esses métodos? É bom ter um link para o artigo, mas preferimos que nossas respostas sejam independentes, especialmente se o link ficar inoperante.
gung - Restabelece Monica
O contexto é a pergunta original: existe multimodalidade? se assim for. e a relevância de um novo método vem do fato de que a caça por impacto de maneira não-paramétrica é um problema de modelagem difícil.
Deep Mukherjee
@gung está pedindo para você expandir sua resposta. Por exemplo, o resultado é de um método explicado em um artigo que não possui versão pública.
Nick Cox
2
Não, quero dizer o que é "LP Nonparametric Mode Identification"? O que é "MaxEnt"? Etc. Em algumas frases, como isso funciona? Por que / quando pode ser preferível a outros métodos? Etc. Estou ciente de que você cria um link para o artigo que os explica, mas seria bom ter algumas frases para apresentá-las aqui, especialmente se o link ficar inoperante, mas mesmo que não para dar aos futuros leitores uma idéia de se eles quer seguir esse método.
gung - Restabelece Monica
2
@DeepMukherjee, você certamente não precisa reescrever o artigo inteiro em sua postagem. Basta adicionar algumas frases dizendo o que é e como funciona.
gung - Restabelece Monica