¿Puedes calcular la potencia de una prueba de Kolmogorov-Smirnov en R?

¿Es posible hacer un análisis de potencia para una prueba de Kolmogorov Smirnov de 2 lados en R?

Estoy probando si dos distribuciones empíricas difieren usando ks.test (), y estoy buscando agregar un análisis de potencia.

No he podido encontrar ningún análisis de potencia incorporado para las pruebas de KS en R. ¿Alguna sugerencia?

Editar : Estas son distribuciones generadas aleatoriamente que se aproximan mucho a mis datos (con tamaños de muestra verdaderos y tasas de disminución estimadas para distribuciones exponenciales)

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)

Estos datos son medidas del tamaño del cuerpo en dos grupos diferentes. Quiero mostrar que los dos grupos tienen esencialmente la misma distribución, pero un colaborador me preguntó si tenía el poder de decir eso en función de los tamaños de muestra. He extraído al azar de una distribución exponencial aquí, pero estos están cerca de los datos reales.

Hasta ahora, he dicho que no hay diferencias significativas en estas distribuciones basadas en la prueba KS de dos lados. También he trazado las dos distribuciones. ¿Cómo puedo demostrar que tengo el poder de hacer tal afirmación, dados los tamaños de muestra y las tasas de desintegración para x e y?

Encontrar poder frente a alternativas de cambio de escala exponencial es razonablemente sencillo.

Sin embargo, no sé si debería usar valores calculados a partir de sus datos para determinar cuál podría ser la potencia. Ese tipo de cálculo de potencia post hoc tiende a dar lugar a conclusiones contrarias a la intuición (y quizás engañosas).

El poder, como el nivel de significancia, es un fenómeno con el que lidias antes del hecho; usaría una comprensión a priori (incluyendo teoría, razonamiento o cualquier estudio previo) para decidir sobre un conjunto razonable de alternativas a considerar y un tamaño de efecto deseable

También puede considerar una variedad de otras alternativas (por ejemplo, podría incorporar el exponencial dentro de una familia gamma para considerar el impacto de casos más o menos sesgados).

Las preguntas habituales que uno podría intentar responder mediante un análisis de poder son:

1) ¿cuál es el poder, para un tamaño de muestra dado, en algún tamaño de efecto o conjunto de tamaños de efecto *?

2) dado un tamaño de muestra y potencia, ¿qué tan grande es detectable un efecto?

3) Dada la potencia deseada para un tamaño de efecto particular, ¿qué tamaño de muestra se requeriría?

* (donde "tamaño del efecto" se pretende genéricamente, y podría ser, por ejemplo, una relación particular de medias, o una diferencia de medias, no necesariamente estandarizada).

Claramente, ya tiene un tamaño de muestra, por lo que no está en el caso (3). Puede considerar razonablemente el caso (2) o el caso (1).

Sugeriría el caso (1) (que también ofrece una forma de tratar el caso (2)).

Para ilustrar un enfoque del caso (1) y ver cómo se relaciona con el caso (2), consideremos un ejemplo específico, con:

• alternativas de cambio de escala

• poblaciones exponenciales

• tamaños de muestra en las dos muestras de 64 y 54

Debido a que los tamaños de muestra son diferentes, tenemos que considerar el caso en que la dispersión relativa en una de las muestras es menor y mayor que 1 (si fueran del mismo tamaño, las consideraciones de simetría hacen posible considerar solo un lado). Sin embargo, debido a que están bastante cerca del mismo tamaño, el efecto es muy pequeño. En cualquier caso, arregle el parámetro para una de las muestras y varíe la otra.

Entonces, lo que uno hace es:

Antemano:

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

Para hacer los cálculos:

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

En R, hice esto:

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

que da la siguiente "curva" de poder

El eje x está en una escala logarítmica, el eje y es la tasa de rechazo.

Es difícil saberlo aquí, pero los puntos negros son ligeramente más altos a la izquierda que a la derecha (es decir, hay una potencia fraccionalmente mayor cuando la muestra más grande tiene la escala más pequeña).

Usando el cdf normal inverso como una transformación de la tasa de rechazo, podemos hacer que la relación entre la tasa de rechazo transformada y log kappa (kappa está sen la gráfica, pero el eje x tiene una escala logarítmica) casi lineal (excepto cerca de 0 ), y el número de simulaciones fue lo suficientemente alto como para que el ruido sea muy bajo; casi podemos ignorarlo para los propósitos actuales.

Entonces podemos usar interpolación lineal. A continuación se muestran los tamaños de efecto aproximados para 50% y 80% de potencia en los tamaños de muestra:

Los tamaños del efecto en el otro lado (el grupo más grande tiene una escala más pequeña) solo se desplazan ligeramente de eso (puede recoger un tamaño de efecto fraccionalmente más pequeño), pero hace poca diferencia, por lo que no trabajaré en el punto.

Entonces, la prueba detectará una diferencia sustancial (de una relación de escalas de 1), pero no una pequeña.

Ahora para algunos comentarios: no creo que las pruebas de hipótesis sean particularmente relevantes para la pregunta de interés subyacente (¿ son bastante similares? ) Y, en consecuencia, estos cálculos de potencia no nos dicen nada directamente relevante para esa pregunta.

Creo que abordas esa pregunta más útil al especificar previamente lo que piensas que "esencialmente lo mismo" realmente significa, operacionalmente. Eso , perseguido racionalmente a una actividad estadística, debería conducir a un análisis significativo de los datos.

Dado que Kolmogorov-Smirnov no es paramétrico, por definición no puede haber un análisis de potencia aplicable. Para tener algún tipo de estimación, debe asumir un modelo de fondo (y, por lo tanto, desviarse del mundo no paramétrico ...) y usarlo para calcular uno de los siguientes: tamaño de muestra, MDE o potencia (es decir, usted arreglar / elegir dos y calcular el tercero).