Valeurs de p égales à 0 dans le test de permutation

J'ai deux jeux de données et je voudrais savoir s'ils sont significativement différents ou non (cela vient de " Deux groupes sont significativement différents? Test à utiliser ").

J'ai décidé d'utiliser un test de permutation, en procédant comme suit dans R:

permutation.test <- function(coding, lncrna) {
    coding <- coding[,1] # dataset1
    lncrna <- lncrna[,1] # dataset2

    ### Under null hyphotesis, both datasets would be the same. So:
    d <- c(coding, lncrna)

    # Observed difference
    diff.observed = mean(coding) - mean(lncrna)
    number_of_permutations = 5000
    diff.random = NULL

    for (i in 1:number_of_permutations) {
        # Sample from the combined dataset
        a.random = sample (d, length(coding), TRUE)
        b.random = sample (d, length(lncrna), TRUE)
        # Null (permuated) difference
        diff.random[i] = mean(b.random) - mean(a.random)
    }

    # P-value is the fraction of how many times the permuted difference is equal or more extreme than the observed difference
    pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations
    pvalue
}

Néanmoins, les valeurs de p ne devraient pas être 0 selon cet article: http://www.statsci.org/smyth/pubs/permp.pdf

Que me recommandez-vous de faire? Est-ce la façon de calculer la valeur de p:

pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations

une bonne façon? Ou vaut-il mieux faire ce qui suit?

pvalue = sum(abs(diff.random) >= abs(diff.observed)) + 1 / number_of_permutations + 1

Discussion

Un test de permutation génère toutes les permutations pertinentes d'un ensemble de données, calcule une statistique de test désignée pour chacune de ces permutations et évalue la statistique de test réelle dans le contexte de la distribution de permutation résultante des statistiques. Une manière courante de l’évaluer est de rendre compte de la proportion de statistiques qui sont (dans un certain sens) "aussi ou plus extrêmes" que les statistiques réelles. Ceci est souvent appelé une «valeur p».

Parce que l'ensemble de données réel est l'une de ces permutations, sa statistique sera nécessairement parmi celles trouvées dans la distribution de permutation. Par conséquent, la valeur de p ne peut jamais être nulle.

Sauf si l'ensemble de données est très petit (moins d'environ 20 à 30 nombres totaux, généralement) ou si la statistique de test a une forme mathématique particulièrement agréable, il n'est pas possible de générer toutes les permutations. (Un exemple où toutes les permutations sont générées apparaît au Test de permutation dans R. ) Par conséquent, les implémentations informatiques des tests de permutation échantillonnent généralement à partir de la distribution de permutation. Ils le font en générant des permutations aléatoires indépendantes et espèrent que les résultats sont un échantillon représentatif de toutes les permutations.

Par conséquent, tous les nombres (tels qu'une "valeur p") dérivés d'un tel échantillon ne sont estimateurs des propriétés de la distribution de permutation. Il est tout à fait possible - et cela arrive souvent lorsque les effets sont importants - que la valeur p estimée soit nulle. Il n'y a rien de mal à cela, mais cela soulève immédiatement la question jusqu'ici négligée de savoir dans quelle mesure la valeur de p estimée pourrait différer de la bonne? Étant donné que la distribution d'échantillonnage d'une proportion (telle qu'une valeur de p estimée) est binomiale, cette incertitude peut être traitée avec un intervalle de confiance binomial .

Architecture

Une mise en œuvre bien conçue suivra de près la discussion à tous égards. Il commencerait par une routine pour calculer la statistique de test, comme celle-ci pour comparer les moyennes de deux groupes:

diff.means <- function(control, treatment) mean(treatment) - mean(control)

Écrivez une autre routine pour générer une permutation aléatoire de l'ensemble de données et appliquez la statistique de test. L'interface avec celle-ci permet à l'appelant de fournir la statistique de test comme argument. Il comparera les premiers méléments d'un tableau (présumé être un groupe de référence) aux éléments restants (le groupe "traitement").

f <- function(..., sample, m, statistic) {
  s <- sample(sample)
  statistic(s[1:m], s[-(1:m)])
}

Le test de permutation est d'abord effectué en trouvant la statistique des données réelles (supposées ici être stockées dans deux tableaux control et treatment), puis en trouvant des statistiques pour de nombreuses permutations aléatoires indépendantes de celles-ci:

z <- stat(control, treatment) # Test statistic for the observed data
sim<- sapply(1:1e4, f, sample=c(control,treatment), m=length(control), statistic=diff.means)

Calculez maintenant l'estimation binomiale de la valeur de p et un intervalle de confiance pour celle-ci. Une méthode utilise le intégrébinconf procédure dans le HMiscpackage:

require(Hmisc)                                    # Exports `binconf`
k <- sum(abs(sim) >= abs(z))                      # Two-tailed test
zapsmall(binconf(k, length(sim), method='exact')) # 95% CI by default

Ce n'est pas une mauvaise idée de comparer le résultat à un autre test, même s'il est connu que ce n'est pas tout à fait applicable: au moins, vous pourriez avoir un ordre de grandeur indiquant où le résultat doit se situer. Dans cet exemple (de comparaison des moyennes), un test t de Student donne généralement un bon résultat de toute façon:

t.test(treatment, control)

Cette architecture est illustrée dans une situation plus complexe, avec des R code de , sur Tester si les variables suivent la même distribution .

Exemple

$10$$0$$20$$1.5$

set.seed(17)
control <- rnorm(10)
treatment <- rnorm(20, 1.5)

Après avoir utilisé le code précédent pour exécuter un test de permutation, j'ai tracé l'échantillon de la distribution de permutation avec une ligne rouge verticale pour marquer la statistique réelle:

h <- hist(c(z, sim), plot=FALSE)
hist(sim, breaks=h$breaks)
abline(v = stat(control, treatment), col="Red")

Le calcul de la limite de confiance binomiale a abouti à

 PointEst Lower        Upper
        0     0 0.0003688199

$0$$0.00037$3.16e-05$0.00037$$0.00037$$0.05$$0.01$$0.001$

commentaires

$k$$N$ $k/N$$(k+1)/(N+1)$$N$

$10$$10^2=100$$0.000005$$1.6$$11.7$parties par million: un peu plus petit que le test t de Student rapporté. Bien que les données aient été générées avec des générateurs de nombres aléatoires normaux, ce qui justifierait l'utilisation du test t de Student, les résultats du test de permutation diffèrent des résultats du test t de Student car les distributions au sein de chaque groupe d'observations ne sont pas parfaitement normales.

$\frac{B}{M}$$\frac{B+1}{M+1}$ ) est un estimateur de valeur p valide (mais conservateur) - il ne conduit pas à un rejet excessif du nul.

(B est le nombre de permutations aléatoires dans lesquelles une statistique supérieure ou égale à celle observée est obtenue et M est le nombre total de permutations aléatoires échantillonnées).

$\frac{B}{M}$