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

15

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
user2886545
la source
(1) La dernière ligne de la question est erronée car elle n'inclut pas les parenthèses nécessaires à l'exécution du calcul prévu. (Il est garanti de produire des résultats supérieurs à , ce qui est impossible pour toute valeur de p.) (2) Vous ne réalisez pas réellement un test de permutation: les deux échantillons et comprendront rarement une partition aléatoire des données mais se chevauchent généralement substantiellement. Au lieu de cela, calculez comme le complément de l' intérieur de l'union de et .1a.randomb.randomb.randoma.randomcodinglncrna
whuber
Parce que la valeur de p est l'ensemble des valeurs au moins aussi extrêmes que celles observées, si l'on évalue la distribution de permutation, la statistique observée est dans les "permutations" comptées. Lors de la randomisation, il est courant de compter la statistique observée parmi les statistiques de permutation considérées (pour des raisons similaires).
Glen_b -Reinstate Monica

Réponses:

15

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

dix0201,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")

Figure

Le calcul de la limite de confiance binomiale a abouti à

 PointEst Lower        Upper
        0     0 0.0003688199

00,000373.16e-050,000370,000370,050,010,001


commentaires

kN k/N(k+1)/(N+1)N

dixdix2=1000,0000051,611,7parties 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.

whuber
la source
L'article de Smyth & Phipson cité ci-dessus montre clairement pourquoi k / N est un mauvais choix pour un estimateur de valeur p. En résumé, pour des niveaux de signification pertinents tels que alpha = 0,05, P ((k / N) <alpha | H0) peut être étonnamment supérieur à alpha. Cela signifie qu'un test de permutation aléatoire utilisant k / N comme estimateur de sa valeur p et 0,05 comme seuil de rejet rejettera l'hypothèse nulle plus de 5% des fois! Une valeur de p nulle est un cas extrême de ce problème - avec un critère alpha = 0, nous nous attendons à ne jamais rejeter le null, mais b / m peut égaler zéro sous le null, conduisant à un faux rejet.
Trisoloriansunscreen
1
@Tal "Un mauvais choix" pour un usage particulier. Ce qui nous distingue en tant que statisticiens des autres, c'est notre compréhension du rôle de la variabilité dans l'analyse des données et la prise de décisions, ainsi que notre capacité à quantifier cette variabilité de manière appropriée. C'est l'approche illustrée (et implicitement préconisée) dans ma réponse ici. Lorsqu'elle est réalisée, il n'y a pas de problème tel que vous le décrivez, car l'utilisateur de la procédure de permutation est amené à comprendre ses limites et ses points forts et aura la liberté d'agir en fonction de ses objectifs.
whuber
13

BMB+1M+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).

BM

Trisoloriansunscreen
la source
1
+1 Ceci est un bon résumé du point principal du document. J'apprécie particulièrement votre attention à la distinction entre une valeur p estimée et une valeur p vraie permutation.
whuber