Teste de randomização / permutação para vetores pareados em R

Eu não sou especialista, então me perdoe se parte da terminologia for um pouco desajeitada. É um prazer fornecer mais informações quando necessário.

Eu tenho dois vetores de 50 valores numéricos emparelhados em R. Eu quero executar um teste de randomização ou permutação bicaudal para determinar se suas diferenças são devidas ao acaso ou não.

Um teste de permutação (também chamado de teste de randomização, teste de re-randomização ou teste exato) é um tipo de teste de significância estatística em que a distribuição da estatística do teste sob a hipótese nula é obtida calculando todos os valores possíveis da estatística do teste sob rearranjos dos rótulos nos pontos de dados observados.

Eu quero fazer esse tipo de teste porque acredito que as distribuições dos valores nos vetores violam as suposições de outros testes, como o teste t (por exemplo, muitos dos valores numéricos no vetor são 0).

A permtestfunção na biblioteca BHH2 quase faz o que eu quero, mas opera em todas as permutações, que levarão muito tempo. Em vez disso, quero estimar o valor p, amostrando um grande número de permutações possíveis. Eu dei uma olhada no pacote de moedas , mas nada parece fazer um teste de permutação com amostras de vetores numéricos emparelhados.$2^{50}$

Alguns pesquisadores me levaram a este e-mail , o que sugere que a razão pela qual não consigo encontrar um pacote para fazer isso é o fato de ser uma linha única no R. Infelizmente, eu não tenho experiência suficiente com o R para poder produzir esse. -forro.

Existe um pacote ou método que executará um teste de permutação emparelhado bicaudal usando apenas uma amostra do espaço de permutação?

Caso contrário, alguém seria capaz de compartilhar um pouco de código R para fazer isso?

Embora eu tenha apontado nos comentários o uso do coinpacote, acho que vale a pena ilustrar que um teste de permutação / randomização é realmente bastante simples, então eu o fiz.

Aqui, escrevo um código R para fazer um teste de randomização para um teste de localização de uma amostra. O teste inverte aleatoriamente os sinais nas diferenças e calcula a média; isso é equivalente a atribuir aleatoriamente cada par de valores aos grupos x e y. O código abaixo pode ser significativamente menor (eu poderia fazê-lo em duas linhas com bastante facilidade, ou mesmo em uma se você não se importasse com código mais lento).

Este código leva alguns segundos na minha máquina:

# assumes the two samples are in 'x' and 'y' and x[i] and y[i] are paired
# set up:
B <- 99999
d <- x-y
m0 <- mean(d)

# perform a one-sample randomization test on d
# for the null hypothesis H0: mu_d = 0   vs H1 mu_d != 0  (i.e. two tailed)
# here the test statistic is the mean
rndmdist <- replicate(B,mean((rbinom(length(d),1,.5)*2-1)*d))

# two tailed p-value:
sum( abs(rndmdist) >= abs(m0))/length(rndmdist)

Essa é a coisa toda.

Observe que rbinom(length(d),1,.5)*2-1)dá um sinal aleatório -1ou 1... ou seja, aleatório; portanto, quando multiplicamos por qualquer conjunto de sinais assinados d, é equivalente a atribuir aleatoriamente +ou -sinais às diferenças absolutas. [Não importa com que distribuição de sinais dvocê comece, agora deles terão sinais aleatórios.]

Aqui, comparo-o com um teste t em alguns dados inventados:

 set.seed(seed=438978)
 z=rnorm(50,10,2)
 x=z-rnorm(50,0,.5)
 y=z+.4+rnorm(50,0,.5)
 t.test(y-x) # gives p = 0.003156

 B <- 99999
 d <- x-y
 m0 <- mean(d)
 rndmdist <- replicate(B,mean((rbinom(length(d),1,.5)*2-1)*d))
 sum( abs(rndmdist) >= abs(m0))/length(rndmdist) 

Quando o teste t é válido, geralmente fornece um valor p muito semelhante ao teste de permutação completamente enumerado, e um valor p simulado como acima (quando o número de simulações é suficientemente grande) convergirá para esse segundo valor p.

No número de repetições usadas acima, um valor p de permutação verdadeiro (ou seja, da enumeração completa) de 0,05 será estimado em 0,001 (ou seja, fornecerá um valor p de randomização entre 0,049 e 0,051) em cerca de 85% do tempo e para 0,002 em 99,5% do tempo.

Aqui está o código para executar um teste de permutação. Eu tenho dados lá, por exemplo. x é a diferença entre os dois vetores.

x <- c(5.1, 9.4, 7.2, 8.1, 8.8, 2.5, 4.2, 6.9, 5.5, 5.3)
m = 5
n = 5
xsum = sum(x)
asum = sum(x[1:m])
bsum = xsum - asum
truediff = asum/m - bsum/n
truediff
abstruediff = abs(truediff)
iter = 100000
difflist <- 1:iter
for(i in 1:iter) {
  s <- sample(x,m) # select a sample of size m
  pasum = sum(s)
  pbsum = sum(x) - sum(s)
  diff  = pasum/m - pbsum/n
  difflist[i] <- diff # add permutation difference to list
}
difflist  <- sort(difflist)
xquantile <- quantile(difflist,probs=c(.005, .01, .025, .05, .95, .975, .99, .995))
xquantile
pdist  <- quantile(difflist, probs=seq(0,1,1/iter))
ntail1 <- length(pdist[difflist <= -abstruediff])
tail1  <- ntail1/iter
tail1  # left-tail probability
ntail2 <- length(pdist[difflist >= abstruediff])
tail2  <- ntail2/iter
tail2  # right-tail probability
twotail = tail1 + tail2
twotail