Como comparar a média de duas amostras cujos dados se encaixam em distribuições exponenciais

10

Eu tenho duas amostras de dados, uma amostra de linha de base e uma amostra de tratamento.

A hipótese é que a amostra de tratamento tenha uma média mais alta que a amostra da linha de base.

Ambas as amostras são de forma exponencial. Como os dados são bastante grandes, só tenho a média e o número de elementos para cada amostra no momento em que executarei o teste.

Como posso testar essa hipótese? Eu acho que é super fácil e já deparei com várias referências ao uso do Teste-F, mas não tenho certeza de como os parâmetros são mapeados.

Jonathan Dobbie
fonte
2
Por que você não tem os dados? Se as amostras são realmente grandes, os testes não paramétricos devem funcionar muito bem, mas parece que você está tentando executar um teste a partir das estatísticas resumidas. Isso está certo?
Mimshot
Os valores de linha de base e tratamento do mesmo conjunto de pacientes ou os dois grupos são independentes?
Michael M
11
@Mimshot, os dados estão fluindo, mas você está certo de que estou tentando executar um teste a partir das estatísticas resumidas. Ele funciona muito bem com um teste Z para dados normais
Jonathan Dobbie
11
Nessas circunstâncias, um teste z aproximado é talvez o melhor que você pode fazer. No entanto, eu me importaria mais com o tamanho do verdadeiro efeito do tratamento, não com a significância estatística. Lembre-se de que, com amostras grandes o suficiente, qualquer pequeno efeito verdadeiro levará a um pequeno valor de p.
Michael M
11
@january - embora, se o tamanho da amostra for grande o suficiente, pelo CLT eles estarão muito próximos da distribuição normal. Sob a hipótese nula, as variações seriam as mesmas (como as médias são); portanto, com um tamanho de amostra grande o suficiente, um teste t deve funcionar bem; não será tão bom quanto você pode fazer com todos os dados, mas ainda assim seria bom. , por exemplo, seria muito bom. n1=n2=100
jbowman

Respostas:

14

Você pode testar a igualdade dos parâmetros médios contra a alternativa de que os parâmetros médios são desiguais com um teste de razão de verossimilhança (teste LR). (No entanto, se os parâmetros médios diferirem e a distribuição for exponencial, será uma mudança de escala, não uma mudança de local.)

Para um teste unilateral (mas apenas assintoticamente no caso bicaudal), acredito que o teste LR é equivalente ao seguinte (para mostrar que, na verdade, é o mesmo que o teste LR para o teste unilateral) caso fosse necessário mostrar que a estatística LR era monotônica em ):x¯/y¯

Digamos que parametrizamos a ésima observação no primeiro exponencial como tendo pdf e a ésima observação no segundo exemplo como tendo pdf (sobre os domínios óbvios para as observações e parâmetros). (Para deixar claro, estamos trabalhando na forma média e não na forma de taxa aqui; isso não afetará o resultado dos cálculos.)1 / μ x exp ( - x i / μ x ) j 1 / μ y exp ( - y j / μ y )i1/μxexp(xi/μx)j1/μyexp(yj/μy)

Como a distribuição de é um caso especial da gama, , a distribuição da soma de 's, é distribuída ; da mesma forma que para a soma dos s, é . Γ ( 1 , μ x ) X S x Γ ( n x , μ x ) Y S y Γ ( n y , μ y )XiΓ(1,μx)XSxΓ(nx,μx)YSyΓ(ny,μy)

Devido ao relacionamento entre distribuições gama e distribuições qui-quadrado, verifica-se que é distribuído . A razão de dois qui-quadrados em seus graus de liberdade é F. Portanto, a razão, .χ 2 2 n x μ y2/μxSxχ2nx2μyμxSx/nxSy/nyF2nx,2ny

Sob a hipótese nula de igualdade de médias, então, e sob a alternativa de dois lados, os valores podem tender a ser menores ou maiores que um valor nulo distribuição, então você precisa de um teste bicaudal.x¯/y¯F2nx,2ny


Simulação para verificar se não cometemos algum erro simples na álgebra:

Aqui simulei 1000 amostras do tamanho 30 para e 20 para partir de uma distribuição exponencial com a mesma média e calculei a estatística da razão de médias acima.YXY

Abaixo está um histograma da distribuição resultante, bem como uma curva que mostra a distribuição que computamos sob o valor nulo:F

exemplo de distribuição simulada da estatística de razão sob o valor nulo


Exemplo, com discussão do cálculo de valores p bicaudais :

Para ilustrar o cálculo, aqui estão duas pequenas amostras de distribuições exponenciais. A amostra X tem 14 observações de uma população com média 10, a amostra Y tem 17 observações de uma população com média 15:

x: 12.173  3.148 33.873  0.160  3.054 11.579 13.491  7.048 48.836 
   16.478  3.323  3.520  7.113  5.358

y:  7.635  1.508 29.987 13.636  8.709 13.132 12.141  5.280 23.447 
   18.687 13.055 47.747  0.334  7.745 26.287 34.390  9.596

As médias da amostra são 12.082 e 16.077, respectivamente. A relação de médias é 0,7515

A área à esquerda é direta, pois está na cauda inferior (calc em R):

 > pf(r,28,34) 
 [1] 0.2210767

Precisamos da probabilidade para a outra cauda. Se a distribuição fosse simétrica no inverso, seria simples fazer isso.

Uma convenção comum com a relação de variâncias do teste F (que é similarmente bicaudal) é simplesmente dobrar o valor p unicaudal (efetivamente o que está acontecendo como aqui ; é também isso que parece ser feito em R, por exemplo ); neste caso, fornece um valor-p de 0,44.

No entanto, se você fizer isso com uma regra de rejeição formal, colocando uma área de em cada cauda, ​​obterá valores críticos conforme descrito aqui . O valor de p é então o maior que levaria à rejeição, o que equivale a adicionar o valor de p unicaudal acima ao valor de p unicaudal na outra cauda para os graus de liberdade trocados. No exemplo acima, que fornece um valor p de 0,43.αα/2α

Glen_b -Reinstate Monica
fonte
Eu estou supondo que este sou apenas eu sendo grosso, mas de onde vem 0,7515?
Jonathan Dobbie
r = média (x) / média (y) = 0,7515 - ou seja, "A proporção de médias"
Glen_b -Reinstala Monica
Ok, incrível. Recebi 0,67, mas isso provavelmente se deve apenas a um erro de entrada de dados.
Jonathan Dobbie
11
Eu fiz a distinção entre as médias da população e a amostra resultante significa mais clara
Glen_b -Reinstate Monica
(+1) Mas, embora seja tangencial, não entendo o último parágrafo. Como dobrar o valor p unicaudal não é equivalente a encontrar o maior , com uma área em cada cauda, ​​que levaria à rejeição? Por que você trocaria os graus de liberdade? ααα2
Scortchi - Restabelece Monica
3

Como um adendo à resposta de @ Glen_b, a taxa de probabilidade é que você pode reorganizar para que . Como existe um mínimo único em , o teste F é realmente o teste da razão de verossimilhança contra alternativas unilaterais à hipótese nula de distribuições idênticas.

nxlognxxi+nylognyyj(nx+ny)lognx+nyxi+yj
nxlog(nxny+1r)+nylog(nynx+r)+nxlognynx+ny+nylognxnx+ny
r=x¯y¯r=1

Para executar o teste de razão de verossimilhança apropriado para uma alternativa de dois lados, você ainda pode usar a distribuição F; você simplesmente precisa encontrar o outro valor da razão da amostra significa para o qual a taxa de probabilidade é igual à da razão observada e . Para este exemplo, , & , fornecendo um valor p geral de (bastante próximo ao obtido pela aproximação do qui-quadrado para a distribuição do dobro da razão de verossimilhança log, ).rELRrobsr E L R = 1,3272 Pr ( R > r E L R ) = 0,2142 0,4352 0,4315Pr(R>rELR)rELR=1.3272Pr(R>rELR)=0.21420.43520.4315

insira a descrição da imagem aqui

Mas dobrar o valor p unicaudal é talvez a maneira mais comum de obter um valor p bicaudal: é equivalente a encontrar o valor da razão da amostra significa para o qual a probabilidade da cauda é igual a e depois encontra . Explicado dessa maneira, pode parecer estar colocando o carro à frente do cavalo para permitir que as probabilidades da cauda definam a extremidade de uma estatística de teste, mas pode ser justificado como sendo efetivamente dois testes de uma cauda (cada um o LRT) com várias comparações correção— & as pessoas geralmente estão interessadas em reivindicar que ou que Pr ( R > r E TrETPPr(R< r o b s )Pr(R> R E T P ) μ x > μ y μ x < μ yPr(R>rETP)Pr(R<robs)Pr(R>rETP)μx>μyμx<μyμ x < μ yμx>μy ou . Também é menos barulhento e, mesmo para amostras de tamanho relativamente pequeno, fornece a mesma resposta que o LRT bicaudal adequado.μx<μy

insira a descrição da imagem aqui

O código R segue:

x <- c(12.173, 3.148, 33.873, 0.160, 3.054, 11.579, 13.491, 7.048, 48.836,
       16.478, 3.323, 3.520, 7.113, 5.358)

y <- c(7.635, 1.508, 29.987, 13.636, 8.709, 13.132, 12.141, 5.280, 23.447, 
       18.687, 13.055, 47.747, 0.334,7.745, 26.287, 34.390, 9.596)

# observed ratio of sample means
r.obs <- mean(x)/mean(y)

# sample sizes
n.x <- length(x)
n.y <- length(y)

# define log likelihood ratio function
calc.llr <- function(r,n.x,n.y){
  n.x * log(n.x/n.y + 1/r) + n.y*log(n.y/n.x + r) + n.x*log(n.y/(n.x+n.y)) + n.y*log(n.x/(n.x+n.y))
}

# observed log likelihood ratio
calc.llr(r.obs,n.x, n.y) -> llr.obs

# p-value in lower tail
pf(r.obs,2*n.x,2*n.y) -> p.lo

# find the other ratio of sample means giving an LLR equal to that observed
uniroot(function(x) calc.llr(x,n.x,n.y)-llr.obs, lower=1.2, upper=1.4, tol=1e-6)$root -> r.hi

#p.value in upper tail
p.hi <- 1-pf(r.hi,2*n.x,2*n.y)

# overall p.value
p.value <- p.lo + p.hi

#approximate p.value
1-pchisq(2*llr.obs, 1)
Scortchi - Restabelecer Monica
fonte