Intervalos de confiança para mediana

9

Eu tenho uma distribuição de amostras com um pequeno número de valores em cada um (menos de ). Eu calculei a mediana para cada amostra, que quero comparar com um modelo e obter a diferença entre o modelo e a mediana de cada amostra. Para ter um resultado consistente, preciso de um erro nessa diferença.10

Isso resulta que encontrar o desvio padrão nesse caso pode ser bastante difícil, pelo menos para um não profissional como eu (veja, por exemplo, aqui ).

Encontrei este site que diz como calcular intervalos de confiança para a mediana, mesmo que não haja referência oficial citada.

Parece-me razoável, mas realmente não posso julgar, então gostaria de saber:

  1. essas fórmulas estão corretas?
  2. Existe uma referência para isso?
  3. E se eu quiser encontrar um IC diferente de ?95%

desde já, obrigado

EDIT: Eu também encontrei este exemplo de bootstrapping para dados não-Gaussianos . Agora, eu não sei muito sobre bootstrapping, mas seria bom ter um endereço sobre sua validade.

Py-ser
fonte
A distribuição exata da amostra de uma mediana da amostra é obtida em stats.stackexchange.com/questions/45124 . (As distribuições assintóticas também são fornecidas na maioria das respostas, mas é improvável que sejam relevantes aqui.) Porém, nenhuma delas é a mesma coisa que um intervalo de confiança ....
whuber
@ whuber, obrigado pelo link, mas não consigo entender a relação. Você poderia ser um pouco mais claro?
Py-ser
Para encontrar um intervalo de confiança (IC) para um parâmetro, usando uma estatística específica, é necessário conhecer a distribuição amostral dessa estatística. Aqui, você procura um IC para a mediana da população (o parâmetro) com base na amostra e pergunta especificamente sobre a mediana da amostra (uma estatística). (O tópico que eu refiro aborda essa última questão.) É crucial conhecer a distribuição exata dessa estatística; daí pode ser derivado um procedimento de intervalo de confiança. Os resultados assintóticos - nos quais sua referência é baseada - correm o risco de serem aproximações ruins para amostras pequenas.
whuber
A estatística é poissoniana. Mas ainda não entendo: a que resultado assintótico você se refere? Essas fórmulas são um caso específico?
Py-ser
1
Acho que você não leu minha resposta nesse tópico , pois fornece um resultado exato para qualquer número de observações: "Essa é uma fórmula exata para a distribuição da mediana para qualquer distribuição contínua".
whuber

Respostas:

14

Sumário

Quando você pode assumir pouco ou nada sobre a verdadeira lei das probabilidades e deduzir pouco sobre ela - como é o caso de pequenas amostras de observações -, um par de estatísticas de ordem adequadamente escolhido constituirá um intervalo de confiança para a mediana. Qual ordem de estatística escolher pode ser facilmente encontrada com uma análise rápida da distribuição Binomial . Existem algumas opções a serem feitas na prática: elas são discutidas e ilustradas no final deste post.( n , 1 / 2 )n(n,1/2)

Aliás, a mesma análise pode ser usada para construir intervalos de confiança para qualquer quantil (do qual a mediana, correspondente a , é um exemplo). A distribuição binomial governa a solução nesse caso.q = 50 % ( n , q )qq=50%(n,q)

Introdução

Lembre-se do que significa um intervalo de confiança (IC). A definição é uma amostra aleatória independente com cada regulada pela mesma distribuição . Supõe-se apenas que é um elemento de um conjunto de possíveis distribuições. Cada um deles tem uma mediana . Para qualquer fixo entre e , um IC de nível é um par de funções (também conhecidas como "estatísticas"), e , de modo queX i F F Ω F 1 / 2 α 0 1 α L LX=(X1,X2,,Xn)XiFFΩF1/2α01αLU

PrF(L(X)F1/2U(X))1α.

O lado direito é a cobertura da CI para a distribuição .F

Além disso: para que isso seja útil, também preferimos que (1) o menor número de coberturas sobre seja o menor possível e (2) a duração esperada do intervalo, , deve tender a ser curto para todos ou "a maioria" .FΩEF(U(X)L(X))FΩ

Análise

Suponha que não assumamos nada sobre . Ω Nesta situação, ainda podemos explorar as estatísticas da ordem . Esses são os valores específicos na amostra classificada. Para simplificar a notação, vamos classificar a amostra de uma vez por todas, para que

X1X2Xn.

O valor é a estatística de ordem da amostra. Como não assumimos nada sobre , não sabemos nada sobre no início, portanto não podemos deduzir muito sobre os intervalos prováveis ​​entre cada e seu vizinho . No entanto, ainda podemos raciocinar quantitativamente sobre os valores individuais: qual é a chance de não exceder a mediana de ? Para descobrir isso, seja uma variável aleatória governada por e permita queXiithΩFXiXi+1XiFYF

πF=PrF(YF1/2)

ser a chance de que não seja superior a mediana de . Então, quando sabemos (desde ) que nossa amostra não ordenada original de valores deve conter pelo menos valores que não excedam .YFXiF1/2X1XiF1/2niF1/2

Este é um problema binomial. Formalmente, se definirmos a variável aleatória como igual a quando e , caso contrário, o anterior mostra que tem uma distribuição de Bernoulli com o parâmetro . Um "sucesso" consiste em observar um valor igual ou inferior à mediana. Portanto, é dada pela probabilidade binomial associada a menos de sucessos:Z1YF1/20ZπFPr(Xi>F1/2)i

Pr(Xi>F1/2)=j=0i1(nj)πFj(1πF)nj.

Você provavelmente notou isso . De fato, para muitas distribuições, os dois valores são iguais: eles diferem apenas quando atribui probabilidade positiva à mediana . Para analisar a diferença, escreva para . Para isso implicaπF1/2FF1/2πF=1/2+εε02(j1)n

πFj(1πF)nj=(1/2+ε)j(1/2ε)nj=(1/2+ε)j[(1/2ε)j(1/2ε)n2j]=(1/4ε2)j(1/2ε)n2j(1/4)j(1/2)n2j=2n.

Consequentemente, quando , podemos nos livrar da dependência da soma de , ao custo de substituir a igualdade por uma desigualdade:2(i1)nF

Pr(Xi>F1/2)2nj=0i1(nj).

Exatamente o mesmo argumento (aplicado pela reversão das estatísticas da ordem) mostra que quando ,2(i+1)n

Pr(Xi<F1/2)2nj=i+1n(nj).

O lado direito reduz a zero sempre que (no primeiro caso) ou (no segundo). Portanto, sempre é possível encontrar índices para os quaisi0inlu

Pr(Xl>F1/2 or Xu<F1/2)=Pr(Xl>F1/2)+Pr(Xu<F1/2)2n(j=0l1(nj)+j=u+1n(nj)).

Solução

Este é o complemento da condição definidora para um intervalo de confiança e, portanto, equivalente a ele:

Pr(XlF1/2Xu)2nj=lu(nj).

Ao selecionar para tornar o lado direito pelo menos , teremos encontrado um procedimento de intervalo de confiança cujo nível é pelo menos .lu1α 1α

Em outras palavras, ao escolher esses índices e , configurando e , o intervalo será um IC para a mediana com cobertura de pelo menos . Você pode calcular sua cobertura real em termos de probabilidades binomiais. Essa cobertura será alcançada para qualquer distribuição que atribua probabilidade zero a (que inclui todas as distribuições contínuas). Será excedido por qualquer que atribua probabilidade diferente de zero a .luL(X)=XlU(X)=Xu[L(X),U(X)]F1/21αFF1/2FF1/2

Discussão

Neste ponto, temos algumas opções. O mais comum é tornar os limites simétricos, definindo razoavelmente próximo de . De fato, estipulando , os limites de confiança podem ser encontrados para qualquer com uma pesquisa rápida ou aplicando a função quantil binomial.un+1lu=n+1ln

Por exemplo, deixe e (para ilustrar um procedimento de CI). Vamos calcular a parte inferior da distribuição binomial cumulativa com os parâmetros e :n=10α=10%1α=90%101/2

> i <- 0:5; names(i) <- i; print(pbinom(i, 10, 1/2), digits=1)
    0     1     2     3     4     5   
0.001 0.011 0.055 0.172 0.377 0.623 

(Esse é um Rcomando e sua resposta.) Como o valor em , igual a , é próximo a , é tentador aceitar e , por então a cobertura será que está próximo da meta de . Se você precisar obter a cobertura desejada, precisará tomar e ou e , ambos com cobertura .25.5%α/2l=3u=10+13=810.0550.055=0.8990%l=2u=8l=3u=910.011.055=0.935

Como verificação, vamos simular muitos conjuntos de dados de qualquer distribuição, calcular esses ICs para os conjuntos de dados e contabilizar a proporção de ICs que cobrem a verdadeira mediana. Este Rexemplo usa uma distribuição Normal:

n <- 10
n.sim <- 1e4
x <- apply(matrix(rnorm(n*n.sim), nrow=n), 2, sort)
covers <- function(x, l, u) mean(x[l, ] <= 0 & x[u, ] >= 0)
c(l3.u8=covers(x,3,8), l2.u8=covers(x,2,8), l3.u9=covers(x,3,9))

A saída é

 l3.u8  l2.u8  l3.u9 
 0.8904 0.9357 0.9319 

As coberturas estão de acordo com os valores teóricos.

Como outro exemplo, vamos desenhar amostras de uma distribuição discreta, como um Poisson:

lambda <- 2
x <- apply(matrix(rpois(n*n.sim, 2), nrow=n), 2, sort)
med <- round(lambda + 1/3 - 0.02/lambda)
c(l3.u8=covers(x,3,8), l2.u8=covers(x,2,8), l3.u9=covers(x,3,9))

 l3.u8  l2.u8  l3.u9 
0.9830 0.9845 0.9964 

Desta vez, as coberturas são muito maiores do que o previsto. O motivo é que há uma chance de que um valor aleatório seja igual à mediana. Isso aumenta muito a chance de o IC cobrir a mediana. Isso não é um problema ou um paradoxo. Por definição, a cobertura deve ser pelo menos não importa qual seja a distribuição - mas é possível (como neste caso) que a cobertura para distribuições específicas seja substancialmente maior que .27%1αF1α

Aí reside a desvantagem: quando você não assume nada sobre , o IC baseado nas estatísticas de pedidos é o único que você pode construir. A cobertura do seu verdadeiro (mas desconhecido) pode ser um pouco maior do que o esperado. Isso significa que seu CI será maior do que se você tinha feito algumas suposições fortes sobre , limitando as possibilidades de .FFΩF

whuber
fonte
Esta resposta se concentra na pergunta nº 3. Quanto às duas primeiras perguntas, (1) ("essas fórmulas estão corretas?"), A resposta não é bem, porque elas usam uma aproximação normal à distribuição binomial; e (2) ("existe uma referência"), a resposta é talvez, mas quem se importa? Uma referência para a análise nesta resposta é Hahn & Meeker, Statistical Interval .
whuber
3

Se você deseja usar métodos numéricos, pode gerar uma estimativa da distribuição amostral de medianas usando o bootstrap. Reamostrar repetidamente sua amostra e calcular muitas medianas. O stdev dessas medianas serve como uma estimativa do stdev da distribuição amostral das medianas. Usei um método semelhante para calcular a incerteza dos resultados dos jogos de xadrez em meu artigo sobre jogos de xadrez que podem ser encontrados aqui https://sonoma.academia.edu/JamalMunshi/papers

Jamal Munshi
fonte
Essa é uma boa ideia. À luz dos comentários à pergunta, o que é necessário é uma análise de sua precisão para pequenos . Além disso, não há motivo para reamostrar repetidamente na prática, porque é fácil obter a distribuição exata de forma fechada. Para um conjunto de dados , a chance de a mediana de uma amostra de bootstrap não exceder (onde ) é a chance de pelo menos metade da os valores de amostra estão no conjunto . Isto é dado por uma distribuição binomial com os parâmetros e . nx1x2xnxxix<xi+1{x1,x2,xi}ni/n
whuber
@ Whuber, desculpe, você quis dizer "isso não é uma boa idéia", certo?
Py-ser
@ Py-ser A idéia subjacente é boa no sentido de que uma versão funcionará, mas a interpretação e a implementação precisam ser aprimoradas.
whuber
Mas toda a nossa discussão anterior foi que você acha que o bootstrapping NÃO é uma boa ideia.
Py-ser