Qual implementação de teste de permutação em R usar em vez de testes t (emparelhados e não emparelhados)?

56

Eu tenho dados de um experimento que analisei usando testes t. A variável dependente é escalonada em intervalos e os dados são desparelhados (ou seja, 2 grupos) ou emparelhados (ou seja, dentro de indivíduos). Por exemplo (dentro dos assuntos):

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

No entanto, como os dados não são normais, um revisor solicitou que usássemos algo diferente do teste t. No entanto, como se pode ver facilmente, os dados não apenas não são normalmente distribuídos, mas também as distribuições não são iguais entre as condições: texto alternativo

Portanto, os testes não paramétricos usuais, o Teste U de Mann-Whitney (não pareado) e o Teste Wilcoxon (emparelhado), não podem ser usados, pois exigem distribuições iguais entre as condições. Por isso, decidi que seria melhor realizar algum teste de reamostragem ou permutação.

Agora, estou procurando uma implementação R de um equivalente baseado em permutação do teste t ou qualquer outro conselho sobre o que fazer com os dados.

Eu sei que existem alguns pacotes R que podem fazer isso por mim (por exemplo, moeda, perm, exactRankTest etc.), mas não sei qual escolher. Portanto, se alguém com alguma experiência com esses testes pudesse me dar um pontapé inicial, isso seria ubercool.

ATUALIZAÇÃO: Seria ideal se você pudesse fornecer um exemplo de como relatar os resultados deste teste.

Henrik
fonte

Respostas:

43

Não deve importar muito, pois a estatística do teste sempre será a diferença de médias (ou algo equivalente). Pequenas diferenças podem advir da implementação dos métodos de Monte-Carlo. Tentando os três pacotes com seus dados com um teste unilateral para duas variáveis ​​independentes:

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
library(coin)                    # for oneway_test(), pvalue()
pvalue(oneway_test(DV ~ IV, alternative="greater", 
                   distribution=approximate(B=9999)))
[1] 0.00330033

library(perm)                    # for permTS()
permTS(DV ~ IV, alternative="greater", method="exact.mc", 
       control=permControl(nmc=10^4-1))$p.value
[1] 0.003

library(exactRankTests)          # for perm.test()
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.003171822

Para verificar o valor p exato com um cálculo manual de todas as permutações, restringirei os dados aos 9 primeiros valores.

x1 <- x1[1:9]
y1 <- y1[1:9]
DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
pvalue(oneway_test(DV ~ IV, alternative="greater", distribution="exact"))
[1] 0.0945907

permTS(DV ~ IV, alternative="greater", exact=TRUE)$p.value
[1] 0.0945907

# perm.test() gives different result due to rounding of input values
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.1029412

# manual exact permutation test
idx  <- seq(along=DV)                 # indices to permute
idxA <- combn(idx, length(x1))        # all possibilities for different groups

# function to calculate difference in group means given index vector for group A
getDiffM <- function(x) { mean(DV[x]) - mean(DV[!(idx %in% x)]) }
resDM    <- apply(idxA, 2, getDiffM)  # difference in means for all permutations
diffM    <- mean(x1) - mean(y1)       # empirical differencen in group means

# p-value: proportion of group means at least as extreme as observed one
(pVal <- sum(resDM >= diffM) / length(resDM))
[1] 0.0945907

coine exactRankTestssão do mesmo autor, mas coinparece ser mais geral e abrangente - também em termos de documentação. exactRankTestsnão é mais desenvolvido ativamente. Portanto, eu escolheria coin(também por causa de funções informativas como support()), a menos que você não goste de lidar com objetos S4.

EDIT: para duas variáveis ​​dependentes, a sintaxe é

id <- factor(rep(1:length(x1), 2))    # factor for participant
pvalue(oneway_test(DV ~ IV | id, alternative="greater",
                   distribution=approximate(B=9999)))
[1] 0.00810081
caracal
fonte
Obrigado pela sua ótima resposta! Mais 2 perguntas: O seu segundo exemplo significa que a moeda realmente fornece toda permutação possível e é um teste exato? Existe algum benefício em não fornecer um teste exato no meu caso?
Henrik
10
(+1) Não é surpresa que o teste t (não emparelhado) produz essencialmente o mesmo valor p, 0,000349. Apesar do que o revisor disse, o teste t é aplicável a esses dados. A razão é que as distribuições amostrais dos meios são aproximadamente normais, mesmo que as distribuições dos dados não sejam. Além disso, como você pode ver pelos resultados, o teste t é realmente mais conservador do que o teste de permutação. (Isto significa que um resultado significativo com o t-teste indica o teste de permutação é susceptível de ser significativo, também.)
whuber
2
@ Henrik Para certas situações (teste escolhido e complexidade numérica), é coinpossível calcular a distribuição exata da permutação (sem realmente passar por todas as permutações, existem algoritmos mais elegantes do que isso). Dada a escolha, a distribuição exata parece preferível, mas a diferença para uma aproximação de Monte Carlo com um número alto de repetições deve ser pequena.
caracal
11
@Caracal Obrigado pelo esclarecimento. Uma pergunta permanece: os dados que apresentei estão emparelhados. Por isso, preciso do equivalente ao teste t emparelhado. A oneway_testfunção é precisa? E se sim, qual é o correto para dados não emparelhados?
Henrik
2
@ Henrik O coinautor me escreveu que oneway_test()não pode calcular a distribuição exata para o caso dependente, você deve seguir a aproximação do MC (somente wilcoxsign_test()é adequado para o teste exato). Eu não sabia disso e preferiria um erro nesse caso, mas o MC deve ser preciso o suficiente com um alto número de repetições.
caracal
29

Alguns comentários estão, acredito, em ordem.

1) Recomendamos que você tente várias exibições visuais de seus dados, porque eles podem capturar coisas que são perdidas por (gráficos como) histogramas, e também recomendo vivamente que plote em eixos lado a lado. Nesse caso, não acredito que os histogramas façam um bom trabalho ao comunicar os principais recursos de seus dados. Por exemplo, dê uma olhada nos boxplots lado a lado:

boxplot(x1, y1, names = c("x1", "y1"))

texto alternativo

Ou mesmo diagramas lado a lado:

stripchart(c(x1,y1) ~ rep(1:2, each = 20), method = "jitter", group.names = c("x1","y1"), xlab = "")

texto alternativo

x1y1x1y1x1y1y1

2) Você não explicou detalhadamente de onde vêm seus dados, nem como eles foram medidos, mas essas informações são muito importantes na hora de selecionar um procedimento estatístico. Suas duas amostras acima são independentes? Existem razões para acreditar que as distribuições marginais das duas amostras devem ser as mesmas (exceto por uma diferença de local, por exemplo)? Quais foram as considerações anteriores ao estudo que o levaram a procurar evidências de uma diferença entre os dois grupos?

zpp

p

5) Na minha opinião, esses dados são um exemplo perfeito (?) De que uma imagem bem escolhida vale 1000 testes de hipóteses. Não precisamos de estatísticas para diferenciar um lápis e um celeiro. A declaração apropriada, a meu ver, para esses dados seria "Esses dados exibem diferenças marcantes em relação à localização, escala e forma". Você pode acompanhar estatísticas descritivas (robustas) de cada uma delas para quantificar as diferenças e explicar o que as diferenças significam no contexto do seu estudo original.

pp

Essa resposta é muito mais longa do que eu pretendia originalmente. Me desculpe por isso.


fonte
Estou curioso para considerar a seguinte abordagem quase visualizada: estimativas de autoinicialização para os momentos dos dois grupos (médias, variações e momentos mais altos, se você desejar) e depois plotar essas estimativas e seus intervalos de confiança, procurando para o grau de sobreposição entre grupos em cada momento. Isso permite falar sobre possíveis diferenças na variedade de características de distribuição. Se os dados estiverem emparelhados, calcule as pontuações da diferença e inicialize os momentos dessa distribuição única. Pensamentos?
Mike Lawrence
2
(+1) Boa análise. Você está perfeitamente certo de que os resultados são óbvios e não é necessário pressionar o ponto com um valor-p. Você pode ser um pouco extremo em sua declaração de (3), porque o teste t não requer dados normalmente distribuídos. Se você estiver preocupado, existem ajustes para a assimetria (por exemplo, a variante de Chen): você pode ver se o valor p para o teste ajustado altera a resposta. Caso contrário, você provavelmente está bem. Nessa situação específica, com esses dados (altamente distorcidos), o teste t funciona bem.
whuber
(+1) Boa captura! E muito bons comentários.
chl
Parece que estamos aceitando a noção de que a distribuição subjacente é "semelhante" à instanciação aleatória. Portanto, não se pode colocar a questão: são ambos do beta (0,25, 0,25) e o teste seria se eles têm o mesmo parâmetro (sem) centralidade. E isso não justificaria o uso de um teste de permutação ou de Wilcoxon?
Dwin
4
5

Meus comentários não são sobre a implementação do teste de permutação, mas sobre as questões mais gerais levantadas por esses dados e sua discussão, em particular o post de G. Jay Kerns.

As duas distribuições são realmente muito parecidas comigo, EXCETO para o grupo de 0s em Y1, que são muito diferentes das outras observações nessa amostra (a menor menor é de cerca de 50 na escala de 0 a 100), bem como todas as do X1. Primeiro, investigaria se havia algo diferente nessas observações.

Segundo, supondo que esses 0s pertençam à análise, dizer que o teste de permutação não é válido porque as distribuições parecem diferir sugere a pergunta. Se o nulo fosse verdadeiro (as distribuições são idênticas), você poderia (com probabilidade razoável) obter distribuições parecendo tão diferentes quanto essas duas? Responder que esse é o objetivo do teste, não é? Talvez neste caso alguns considerem a resposta óbvia sem executar o teste, mas com essas distribuições pequenas e peculiares, acho que não.

Andrew Taylor
fonte
Parece que este deve ser um ou mais comentários, não uma resposta. Se você clicar no pequeno "adicionar comentário" cinza, poderá colocar seus pensamentos na conversa abaixo da pergunta ou de uma resposta específica, onde eles pertencem. Você faz pontos substantivos aqui, mas não está claro que este não seja o local apropriado para eles.
gung - Restabelece Monica
11
@gung É preciso um pouco de reputação para poder postar um comentário ;-).
whuber
4
2(1/2)7.01
4

Como essa pergunta apareceu novamente, posso adicionar outra resposta inspirada em uma postagem recente do blog via R-Bloggers de Robert Kabacoff, autor do Quick-R e R in Action usando o lmPermpacote.

No entanto, esse método produz resultados nitidamente contrastantes (e muito instáveis) aos produzidos pelo coinpacote na resposta de @caracakl (o valor p da análise dentro dos sujeitos é 0.008). A análise também toma a preparação dos dados da resposta de @ caracal:

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
id <- factor(rep(1:length(x1), 2)) 

library(lmPerm)

summary(aovp( DV ~ IV + Error(id)))

produz:

> summary(aovp( DV ~ IV + Error(id)))
[1] "Settings:  unique SS "

Error: id
Component 1 :
          Df R Sum Sq R Mean Sq
Residuals 19    15946       839


Error: Within
Component 1 :
          Df R Sum Sq R Mean Sq Iter Pr(Prob)  
IV         1     7924      7924 1004    0.091 .
Residuals 19    21124      1112                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Se você executar isso várias vezes, os valores p saltam entre ~ 0,05 e ~ .1.

Embora seja uma resposta para a pergunta, permita-me fazer uma pergunta no final (eu posso passar para uma nova pergunta, se desejar):
Qualquer idéia de por que essa análise é tão instável e produz valores p tão divergentes para a análise de moedas? Fiz algo de errado?

Henrik
fonte
2
Pode ser melhor fazer isso como uma pergunta separada, se realmente é uma pergunta que você deseja responder, em vez de outra solução possível que você deseja listar com o restante. Percebo que você especifica um estrato de erro, mas @caracal não; esse seria meu primeiro palpite sobre a diferença entre b / t e este. Além disso, ao simular, os valores geralmente saltam; para reprodutibilidade, você especifica a semente, por exemplo set.seed(1); para maior precisão na estimativa de MC, você aumenta o número de iterações; Não tenho certeza se uma dessas é a resposta 'certa' para sua pergunta, mas provavelmente é relevante.
gung - Restabelece Monica
2
Novamente, sugiro comparar os resultados do MC com os cálculos manuais usando o teste de permutação total (re-randomização). Veja o código do seu exemplo para obter uma comparação de oneway_anova()(sempre perto do resultado correto) e aovp()(geralmente longe do resultado correto). Não sei por que aovp()dá resultados bastante variados, mas pelo menos nesse caso aqui são implausíveis. @gung a última chamada oneway_test(DV ~ IV | id, ...)na minha resposta original especificou os estratos de erro para o caso dependente (sintaxe diferente da usada por aov()).
Caracal
@caracal, você está certo. Não olhei para o último bloco de código após a edição. Eu estava olhando para o bloco de código superior - desleixado da minha parte.
gung - Restabelece Monica
Eu realmente não preciso da resposta. É apenas mais uma possibilidade que vale a pena mencionar aqui. Infelizmente, está longe de outros resultados, os quais também vale a pena notar.
Henrik
11
@Henrik executa o Aovp com maxExact = 1000. Se demorar muito, defina iter = 1000000 e Ca = 0.001. O cálculo termina quando o erro padrão estimado de p é menor que Ca * p. (Os valores mais baixos dão resultados mais estáveis.)
xmjx
1

Essas pontuações são proporções? Nesse caso, você certamente não deveria estar usando um teste paramétrico gaussiano e, embora possa avançar com uma abordagem não-paramétrica, como um teste de permutação ou auto-inicialização dos meios, sugiro que você obtenha mais poder estatístico empregando uma abordagem paramétrica não gaussiana adequada. Especificamente, sempre que você puder calcular uma medida de proporção dentro de uma unidade de interesse (por exemplo, participante de um experimento), poderá e provavelmente deverá usar um modelo de efeitos mistos que especifique observações com erro distribuído binomialmente. Veja Dixon 2004 .

Mike Lawrence
fonte
As pontuações não são proporções, mas estimativas dos participantes em uma escala de 0 a 100 (os dados apresentados são médias de estimativas em vários itens dessa escala).
Henrik
Então os não paramétricos pareceriam o caminho tradicional a seguir. Dito isso, perguntei-me se esses dados de escala poderiam ser úteis para derivar de um processo binomial e, assim, analisados ​​como tal. Ou seja, você diz que cada pontuação é a média de vários itens, e digamos que cada item é uma escala de 10 pontos; nesse caso, eu representaria uma resposta de, digamos, "8" como uma série de tentativas, 8 de que possuem o valor 1 e dois com o valor 0, todos rotulados com o mesmo rótulo em uma variável "item". Com esses dados expandidos / binomiais, você pode calcular o modelo de efeitos mistos binomiais.
Mike Lawrence
Seguindo o meu comentário anterior, devo observar que nos dados expandidos / binomiais, você pode modelar a variável "item" como um efeito fixo ou aleatório. Eu acho que preferiria modelá-lo como um efeito fixo, porque presumivelmente você pode estar interessado não apenas em contabilizar, mas também em avaliar as diferenças de itens e qualquer possível interação entre o item e outras variáveis ​​preditivas.
Mike Lawrence
0

Apenas adicionando outra abordagem, ezPermdo ezpacote:

> # preparing the data
> DV <- c(x1, y1)
> IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
> id <- factor(rep(1:length(x1), 2))
> df <- data.frame(id=id,DV=DV,IV=IV)
>
> library(ez)
> ezPerm( data = df, dv = DV, wid = id, within = IV, perms = 1000)
|=========================|100%              Completed after 17 s 
  Effect     p p<.05
1     IV 0.016     *

Isso parece ser consistente com o oneway_testdo coinpacote:

> library(coin)
> pvalue(oneway_test(DV ~ IV | id,  distribution=approximate(B=999999)))
[1] 0.01608002
99 percent confidence interval:
 0.01575782 0.01640682

No entanto, observe que este não é o mesmo exemplo fornecido pelo @caracal . Em seu exemplo, ele inclui alternative="greater", portanto, a diferença de valores de p ~0.008vs ~0.016.

O aovppacote sugerido em uma das respostas produzem desconfiado reduzir os valores de p, e corre de forma suspeita rápido mesmo quando eu tento altos valores para a Iter, Cae maxIterargumentos:

library(lmPerm)
summary(aovp(DV ~ IV + Error(id/IV), data=df,  maxIter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Iter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Ca = 0.00000000001))

Dito isto, os argumentos parecem reduzir um pouco as variações dos valores de p ~.03e ~.1(obtive um intervalo maior do que o relatado por @Henrik), para 0.03e 0.07.

toto_tico
fonte