Pouvez-vous calculer la puissance d'un test de Kolmogorov-Smirnov dans R?

Est-il possible de faire une analyse de puissance pour un test de Kolmogorov Smirnov bilatéral en R?

Je teste si deux distributions empiriques diffèrent en utilisant ks.test () et cherche à ajouter une analyse de puissance.

Je n'ai pas pu trouver d'analyses de puissance intégrées pour les tests KS dans R. Des suggestions?

Edit : Ce sont des distributions générées de manière aléatoire qui se rapprochent étroitement de mes données (avec des tailles d'échantillon réelles et des taux de décroissance estimés pour des distributions exponentielles)

set.seed(100)
x <- rexp(64, rate=0.34)
y <- rexp(54,rate=0.37)

#K-S test: Do x and y come from same distribution?
ks.test(x,y)

Ces données sont des mesures de la taille corporelle dans deux groupes différents. Je veux montrer que les deux groupes ont essentiellement la même distribution, mais un collaborateur m'a demandé si j'avais le pouvoir de le dire en fonction de la taille des échantillons. J'ai tiré au hasard d'une distribution exponentielle ici, mais celles-ci sont proches des données réelles.

Jusqu'à présent, j'ai dit qu'il n'y avait pas de différence significative dans ces distributions sur la base du test KS bilatéral. J'ai également tracé les deux distributions. Comment puis-je montrer que j'ai le pouvoir de faire une telle déclaration, compte tenu de la taille des échantillons et des taux de décroissance pour x et y?

Trouver du pouvoir contre des alternatives exponentielles de changement d'échelle est assez simple.

Cependant, je ne sais pas si vous devez utiliser des valeurs calculées à partir de vos données pour déterminer quelle pourrait être la puissance. Ce type de calcul de puissance post hoc a tendance à aboutir à des conclusions contre-intuitives (et peut-être trompeuses).

Le pouvoir, comme le niveau de signification, est un phénomène que vous traitez avant le fait; vous utiliseriez une compréhension a priori (y compris la théorie, le raisonnement ou toute étude antérieure) pour décider d'un ensemble raisonnable d'alternatives à considérer et d'une taille d'effet souhaitable

Vous pouvez également envisager une variété d'autres alternatives (par exemple, vous pouvez intégrer l'exponentielle à l'intérieur d'une famille gamma pour prendre en compte l'impact de cas plus ou moins asymétriques).

Les questions habituelles auxquelles on pourrait essayer de répondre par une analyse de puissance sont:

1) quelle est la puissance, pour une taille d'échantillon donnée, pour une certaine taille d'effet ou un ensemble de tailles d'effet *?

2) étant donné la taille et la puissance de l'échantillon, quelle est la taille d'un effet détectable?

3) Étant donné une puissance souhaitée pour une taille d'effet particulière, quelle taille d'échantillon serait requise?

* (où ici la «taille d'effet» est conçue de manière générique et peut être, par exemple, un rapport particulier de moyennes, ou une différence de moyennes, pas nécessairement standardisée).

De toute évidence, vous avez déjà une taille d'échantillon, vous n'êtes donc pas dans le cas (3). Vous pouvez raisonnablement envisager le cas (2) ou le cas (1).

Je suggère le cas (1) (qui donne également un moyen de traiter le cas (2)).

Pour illustrer une approche du cas (1) et voir comment elle se rapporte au cas (2), considérons un exemple spécifique, avec:

• alternatives de changement d'échelle

• populations exponentielles

• tailles d'échantillon dans les deux échantillons de 64 et 54

Parce que les tailles d'échantillon sont différentes, nous devons considérer le cas où l'écart relatif dans l'un des échantillons est à la fois plus petit et plus grand que 1 (s'ils étaient de la même taille, des considérations de symétrie permettent de ne considérer qu'un seul côté). Cependant, comme ils sont assez proches de la même taille, l'effet est très faible. Dans tous les cas, fixez le paramètre pour l'un des échantillons et modifiez l'autre.

Donc ce que l'on fait c'est:

Préalablement:

choose a set of scale multipliers representing different alternatives
select an nsim (say 1000)
set mu1=1

Pour effectuer les calculs:

for each possible scale multiplier, kappa 
  repeat nsim times
    generate a sample of size n1 from Exp(mu1) and n2 from Exp(kappa*mu1)
    perform the test
  compute the rejection rate across nsim tests at this kappa

Dans R, j'ai fait ceci:

alpha = 0.05
n1 = 54
n2 = 64
nsim = 10000
s = c(1.1,1.2,1.5,2,2.5,3) # set up grid for kappa
s = c(1/rev(s),1,s)        #  also below and at 1
rr = array(NA,length(s))   # to hold rejection rates

for(i in seq_along(s)) rr[i]=mean(replicate(nsim,
                                    ks.test(rexp(n1,1),rexp(n2,s[i]))$p.value)<alpha
                                 )

plot(rr~s,log="x",ylim=c(0,1),type="n") #set up plot
points(rr~rev(s),col=3) # plot the reversed case to show the (tiny) asymmetry+noise
points(rr~s,col=1) # plot the "real" case last 
abline(h=alpha,col=8,lty=2) # draw in alpha

ce qui donne la "courbe" de puissance suivante

L'axe des x est sur une échelle logarithmique, l'axe des y est le taux de rejet.

C'est difficile à dire ici, mais les points noirs sont légèrement plus élevés sur la gauche que sur la droite (c'est-à-dire qu'il y a une fraction de plus de puissance lorsque le plus grand échantillon a la plus petite échelle).

En utilisant le cdf normal inverse comme transformation du taux de rejet, nous pouvons faire la relation entre le taux de rejet transformé et log kappa (kappa est sdans le graphique, mais l'axe des x est à l'échelle logarithmique) très presque linéaire (sauf près de 0 ), et le nombre de simulations était suffisamment élevé pour que le bruit soit très faible - nous pouvons à peu près l'ignorer aux fins actuelles.

Nous pouvons donc simplement utiliser l'interpolation linéaire. Vous trouverez ci-dessous des tailles d'effet approximatives pour une puissance de 50% et 80% à vos tailles d'échantillon:

Les tailles d'effet de l'autre côté (un groupe plus grand a une échelle plus petite) ne sont que légèrement décalées par rapport à cela (peuvent prendre une taille d'effet légèrement plus petite), mais cela ne fait pas beaucoup de différence, donc je ne travaillerai pas sur le point.

Ainsi, le test détectera une différence substantielle (à partir d'un rapport d'échelles de 1), mais pas une petite.

Maintenant, pour quelques commentaires: je ne pense pas que les tests d'hypothèse soient particulièrement pertinents pour la question d'intérêt sous-jacente ( sont-ils assez similaires? ), Et par conséquent ces calculs de puissance ne nous disent rien de directement lié à cette question.

Je pense que vous répondez à cette question plus utile en préspécifiant ce que vous pensez «essentiellement le même» signifie réellement, sur le plan opérationnel. Cela - poursuivi rationnellement en une activité statistique - devrait conduire à une analyse significative des données.

Étant donné que Kolmogorov-Smirnov n'est pas paramétrique, par définition, il ne peut y avoir d'analyse de puissance applicable. Pour avoir une sorte d'estimation, vous devez supposer un modèle d'arrière-plan (et donc vous détourner du monde non paramétrique ...) et l'utiliser pour calculer l'un des éléments suivants: taille de l'échantillon, MDE ou puissance (c'est-à-dire, vous fixer / choisir deux et calculer le troisième).