Testez si les variables suivent la même distribution

Si vous voulez tester si deux variables suivent la même distribution, serait-ce un bon test de simplement trier les deux variables, puis de vérifier leur corrélation? S'il est élevé (au moins 0,9?), Les variables proviennent probablement de la même distribution.

Par distribution, je veux dire "normal", "chi carré", "gamma" etc.

Voyons si c'est un bon test ou non. Il y a beaucoup plus que simplement prétendre que c'est mauvais ou montrer dans un cas que cela ne fonctionne pas bien. La plupart des tests fonctionnent mal dans certaines circonstances, si souvent nous sommes confrontés à l'identification des circonstances dans lesquelles tout test proposé pourrait éventuellement être un bon choix.

Description du test

Comme tout test d'hypothèse, celui-ci consiste en (a) une hypothèse nulle et alternative et (b) une statistique de test (le coefficient de corrélation) destinée à discriminer les hypothèses.

L'hypothèse nulle est que les deux variables proviennent de la même distribution. Pour être précis, nommons les variables $X$ et $Y$ et supposons que nous avons observé $n_x$ instances de $X$ , appelées $x_i = (x_1, x_2, \ldots, x_{n_x})$ , et $n_y$ instances de $Y$ , appelées $y_i$ . L'hypothèse nulle est que toutes les instances de $X$ et $Y$ sont indépendantes et identiquement distribuées (iid).

Prenons comme hypothèse alternative que (a) toutes les instances de sont iid selon une distribution sous-jacente F X et (b) toutes les instances de Y sont iid selon une distribution sous-jacente F Y mais (c) F X diffère de F Oui . (Ainsi, nous ne chercherons pas les corrélations entre les x i , les corrélations entre les y i , les corrélations entre les x i et y j , ni les différences de distribution entre les x ou $X$$F_X$$Y$$F_Y$$F_X$$F_Y$$x_i$$y_i$$x_i$$y_j$$x$$y$séparément: on suppose que ce n'est pas plausible.)

La statistique de test proposée suppose que (appelez cette valeur commune n ) et calcule le coefficient de corrélation de ( x [ i ] , y [ i ] ) (où, comme d'habitude, [ i ] désigne le i ème plus petit des données). Appelons cela t ( x , y ) .$n_x = n_y$$n$$(x_{[i]}, y_{[i]})$$[i]$$i^\text{th}$$t(x,y)$

Tests de permutation

Dans cette situation - quelle que soit la statistique proposée - nous pouvons toujours effectuer un test de permutation. Sous l'hypothèse nulle, la probabilité des données ( ( x 1 , x 2 , … , x n ) , ( y 1 , y 2 , … , y n ) ) est la même que la probabilité d'une permutation des 2 n valeurs de données. En d'autres termes, l'affectation de la moitié des données à X et de l'autre moitié à $t$$\left((x_1, x_2, \ldots, x_n), (y_1, y_2, \ldots, y_n)\right)$$2n$$X$$Y$est une pure coïncidence aléatoire. Ceci est une simple conséquence directe des hypothèses iid et l'hypothèse nulle que .$F_X=F_Y$

Par conséquent, la distribution d'échantillonnage de , conditionnée aux observations x i et y i , est la distribution de toutes les valeurs de t atteintes pour tous ( 2 n ) ! permutations des données. Nous nous intéressons à cela parce que pour toute taille de test prévue donnée α , telle que α = 0,05 (correspondant à une confiance de 95 %), nous construirons une région critique bilatérale à partir de la distribution d'échantillonnage de t : elle se compose de la plus extrême$t(x,y)$$x_i$$y_i$$t$$(2n)!$$\alpha$$\alpha = .05$$95$$t$ % des valeurs possibles de t (sur le côté haut, car une corrélation élevée est compatible avec des distributions similaires et une faible corrélation ne l'est pas). C'est ainsi que nous déterminons la taille du coefficient de corrélation afin de décider que les données proviennent de différentes distributions.$100\alpha$$t$

Simulation de la distribution d'échantillonnage nulle

Parce que (ou, si vous voulez, ( 2 $(2n)!$, qui compte le nombre de façons de diviser les2ndonnées en deux morceaux de taillen) devient grand même pour les petitsn, il n'est pas possible de calculer exactement la distribution d'échantillonnage, nous l'échantillons donc à l'aide d'une simulation. (Par exemple, lorsquen=16, ( 2$\binom{2n}{n}/2$$2n$$n$$n$$n=16$et(2n)! ≈2,63×1035.) Un millier d'échantillons suffit souvent (et certainement pour les explorations que nous allons entreprendre).$\binom{2n}{n}/2 = 300\ 540\ 195$$(2n)! \approx 2.63\times 10^{35}$

Il y a deux choses que nous devons découvrir: premièrement, à quoi ressemble la distribution d'échantillonnage sous l'hypothèse nulle. Deuxièmement, dans quelle mesure ce test établit-il une distinction entre les différentes distributions?

Il y a une complication: la distribution d'échantillonnage dépend de la nature des données. Tout ce que nous pouvons faire est de regarder des données réalistes, créées pour émuler tout ce que nous sommes intéressés à étudier, et espérons que ce que nous apprenons des simulations s'appliquera à notre propre situation.

la mise en oeuvre

Pour illustrer, j'ai réalisé ce travail en R. Il tombe naturellement en trois morceaux.

1. Une fonction pour calculer la statistique de test . Parce que je veux être un peu plus général, ma version gère des ensembles de données de tailles différentes ( n x ≠ n y ) en interpolant linéairement parmi les valeurs de l'ensemble de données plus grand (trié) pour créer des correspondances avec l'ensemble de données plus petit (trié). Parce que cela est déjà fait par la fonction , je prends juste ses résultats:$t(x,y)$$n_x \ne n_y$Rqqplot

   test.statistic <- function(x, y) {
     transform <- function(z) -log(1-z^2)/2
     fit <- qqplot(x,y, plot.it=FALSE)
     transform(cor(fit$x, fit$y))
   }

   Une petite torsion - inutile mais utile pour la visualisation - ré-exprime le coefficient de corrélation d'une manière qui rendra la distribution de la statistique nulle approximativement symétrique. Voilà ce qui transformse passe.

2. La simulation de la distribution d'échantillonnage. Pour l'entrée, cette fonction accepte le nombre d'itérations n.iterainsi que les deux ensembles de données dans les tableaux xet y. Il génère un tableau de n.itervaleurs de la statistique de test. Son fonctionnement interne doit être transparent, même pour un non- Rutilisateur:

   permutation.test <- function(n.iter, x, y) {
     z <- c(x,y)
     n.x <- length(x)
     n.y <- length(y)
     n <- length(z)
     k <- min(n.x, n.y)
     divide <- function() {
       i <- sample.int(n, size=k)
       test.statistic(z[i], z[-i])
     }
     replicate(n.iter, divide())
   }

3. Bien que ce soit tout ce dont nous avons besoin pour effectuer le test, afin de l'étudier, nous voudrons répéter le test plusieurs fois. Ainsi, nous effectuons le test une fois et enveloppons ce code dans une troisième couche fonctionnelle, généralement nommée fici, que nous pouvons appeler à plusieurs reprises. Pour le rendre suffisamment général pour une étude large, il accepte en entrée la taille des ensembles de données à simuler ( n.xet n.y), le nombre d'itérations pour chaque test de permutation ( n.iter), une référence à la fonction testpour calculer la statistique de test (vous verrez momentanément pourquoi nous pourrions ne pas vouloir coder en dur ceci), et deux fonctions pour générer des valeurs aléatoires iid, une pour ( ) et une pour Y ( ). Une option$X$dist.x$Y$dist.yplot.it est utile pour voir ce qui se passe.

   f <- function(n.x, n.y, n.iter, test=test.statistic, dist.x=runif, dist.y=runif, 
       plot.it=FALSE) {
     x <- dist.x(n.x)
     y <- dist.y(n.y)
     if(plot.it) qqplot(x,y)
   
     t0 <- test(x,y)
     sim <- permutation.test(n.iter, x, y)
     p <- mean(sim > t0) + mean(sim==t0)/2
     if(plot.it) {
       hist(sim, xlim=c(min(t0, min(sim)), max(t0, max(sim))), 
            main="Permutation distribution")
       abline(v=t0, col="Red", lwd=2)
     }
     return(p)
   }

   La sortie est une "valeur p" simulée: la proportion de simulations produisant une statistique qui semble plus extrême que celle réellement calculée pour les données.

Les parties (2) et (3) sont extrêmement générales: vous pouvez effectuer une étude comme celle-ci pour un test différent en remplaçant simplement par test.statisticun autre calcul. Nous le faisons ci-dessous.

Premiers résultats

Par défaut, notre code compare les données tirées de deux distributions uniformes. Je le laisse faire (pour , qui sont des ensembles de données assez petits et présentent donc un cas de test modérément difficile), puis je le répète pour une comparaison uniforme-normale et une comparaison uniforme-exponentielle. (Les distributions uniformes ne sont pas faciles à distinguer des distributions normales à moins que vous ayez un peu plus de 16 valeurs, mais les distributions exponentielles - ayant une asymétrie élevée et une longue queue droite - se distinguent généralement facilement des distributions uniformes.)$n.x = n.y = 16$$16$

set.seed(17)             # Makes the results reproducible
n.per.rep <- 1000        # Number of iterations to compute each p-value
n.reps <- 1000           # Number of times to call `f`
n.x <- 16; n.y <- 16     # Dataset sizes

par(mfcol=c(2,3))        # Lay results out in three columns
null <- replicate(n.reps, f(n.x, n.y, n.per.rep))
hist(null, breaks=20)
plot(null)

normal <- replicate(n.reps, f(n.x, n.y, n.per.rep, dist.y=rnorm))
hist(normal, breaks=20)
plot(normal)

exponential <- replicate(n.reps, f(n.x, n.y, n.per.rep, dist.y=function(n) rgamma(n, 1)))
hist(exponential, breaks=20)
plot(exponential)

$X$$Y$$X$$Y$

$16$$x_i$$16$$y_i$f$0.05$$11$$16$des valeurs indépendantes de chacun. C'est une puissance assez faible. Mais c'est peut-être inévitable, alors continuons.

Les graphiques de droite testent de la même manière une distribution uniforme contre une distribution exponentielle. Ce résultat est bizarre. Ce test tend, le plus souvent, à conclure que les données uniformes et les données exponentielles se ressemblent. Il semble "penser" que les variables uniformes et exponentielles sont plus similaires que deux variables uniformes! Que se passe t-il ici?

Le problème est que les données d'une distribution exponentielle auront tendance à avoir quelques valeurs extrêmement élevées. Lorsque vous faites un nuage de points de ceux-ci par rapport à des valeurs uniformément distribuées, il y aura alors quelques points loin en haut à droite de tout le reste. Cela correspond à un coefficient de corrélation très élevé. Ainsi, chaque fois que l'une ou l'autre des distributions génère quelques valeurs extrêmes, le coefficient de corrélation est un choix terrible pour mesurer la différence entre les distributions. Cela conduit à un autre problème encore pire: à mesure que la taille des ensembles de données augmente, les chances d'obtenir quelques observations extrêmes augmentent. Ainsi, nous pouvons nous attendre à ce que ce test soit de plus en plus mauvais à mesure que la quantité de données augmente. Comme c'est horrible ....

Un meilleur test

$y=x$

Voici une Rimplémentation:

test.statistic <- function(x, y) {
  ks.test(x,y)$statistic
}

C'est vrai: il est intégré au logiciel, nous n'avons donc qu'à l'appeler. Mais attendez! Si vous lisez attentivement le manuel, vous apprendrez que (a) le test fournit une valeur de p mais (b) que la valeur de p est (grossièrement) incorrecte lorsque les deux xet ysont des ensembles de données. Il est destiné à être utilisé lorsque vous pensez savoir exactement de quelle distribution xproviennent les données et que vous voulez voir si c'est vrai. Ainsi, le test ne tient pas correctement compte de l'incertitude quant à la distribution des données y.

Aucun problème! Le cadre de test de permutation est toujours aussi valide. En effectuant le changement précédent sur test.statistic, tout ce que nous avons à faire est de relancer l'étude précédente, sans changement. Voici les résultats.

$p=0.20$

$70$$0.05$$11$

$30$$\alpha=5$$50$$\alpha=10$$0.10$

Conclusions

Ainsi, les problèmes avec le test de corrélation ne sont pas dus à une difficulté inhérente à ce paramètre. Non seulement le test de corrélation fonctionne très mal, mais il est mauvais par rapport à un test largement connu et disponible. (Je suppose qu'elle est inadmissible, ce qui signifie qu'elle fonctionnera toujours moins bien, en moyenne, que la version à permutation du test KS, ce qui implique qu'il n'y a aucune raison de l'utiliser.)

Non, la corrélation n'est pas un bon test de cela.

x <- 1:100 #Uniform
y <- sort(rnorm(100)) #Normal
cor(x,y) #.98

Je ne connais pas de bon test qui compare si, par exemple, deux distributions sont à la fois normales, mais éventuellement avec une moyenne et une sd différentes Indirectement, vous pouvez tester la normalité de chacune, séparément, et si les deux semblaient normales, devinez qu'elles étaient toutes les deux.

S'il y a un nombre suffisamment grand de variables, cela peut montrer plus de corrélation avec les valeurs classées par taille. Cependant, elle ne semble pas être une méthode particulièrement utile, notamment parce qu'elle fournit peu de moyens d'estimer la confiance qu'ils pourraient utiliser le même modèle.

Un problème que vous êtes susceptible de rencontrer est lorsque vous avez des modèles avec une moyenne et une asymétrie similaires, mais une différence de kurtosis, car un nombre modéré de mesures peut convenir suffisamment pour avoir l'air assez bien corrélé.

Il semble plus raisonnable de modéliser les deux variables par rapport à des distributions différentes pour voir laquelle est la plus probable pour chacune et comparer les résultats.

Il peut être utile de normaliser les deux valeurs, de trier et de tracer chacune - cela vous permettra de voir comment les ajustements se comparent - et vous pouvez également tracer un modèle possible pour les deux, qui serait lié à ce que vous avez suggéré, mais plutôt que en attendant une réponse concrète, juste une idée visuelle sur la proximité des distributions.