Você pode calcular a potência de um teste de Kolmogorov-Smirnov em R?

É possível fazer uma análise de potência para um teste Kolmogorov Smirnov de dois lados em R?

Estou testando se duas distribuições empíricas diferem usando ks.test () e pretendo adicionar uma análise de potência.

Não consegui encontrar nenhuma análise de potência integrada para testes KS em R. Alguma sugestão?

Editar : são distribuições geradas aleatoriamente que aproximam meus dados (com tamanhos de amostra verdadeiros e taxas de decaimento estimadas para distribuições exponenciais)

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)

Esses dados são medidas do tamanho do corpo em dois grupos diferentes. Quero mostrar que os dois grupos têm essencialmente a mesma distribuição, mas um colaborador me perguntou se eu tinha o poder de dizer isso com base no tamanho da amostra. Eu desenhei aleatoriamente a partir de uma distribuição exponencial aqui, mas estas são próximas aos dados reais.

Até agora, eu disse que não há diferença significativa nessas distribuições com base no teste KS nos dois lados. Eu também plotei as duas distribuições. Como posso mostrar que tenho o poder de fazer tal afirmação, dados os tamanhos das amostras e as taxas de decaimento para x e y?

Encontrar poder contra alternativas exponenciais de mudança de escala é razoavelmente simples.

No entanto, não sei se você deve usar valores calculados a partir de seus dados para descobrir qual poderia ter sido o poder. Esse tipo de cálculo de poder post hoc tende a resultar em conclusões contra-intuitivas (e talvez enganosas).

O poder, como o nível de significância, é um fenômeno com o qual você lida antes do fato; você usaria um entendimento a priori (incluindo teoria, raciocínio ou estudos anteriores) para decidir sobre um conjunto razoável de alternativas a considerar e um tamanho de efeito desejável

Você também pode considerar uma variedade de outras alternativas (por exemplo, você pode incorporar o exponencial dentro de uma família gama para considerar o impacto de casos mais ou menos assimétricos).

As perguntas comuns que se pode tentar responder por uma análise de poder são:

1) qual é o poder, para um determinado tamanho de amostra, em algum tamanho ou conjunto de tamanhos de efeito *?

2) dado um tamanho e poder de amostra, qual o tamanho de um efeito é detectável?

3) Dada a potência desejada para um tamanho de efeito específico, qual tamanho de amostra seria necessário?

* (onde aqui 'tamanho do efeito' se destina genericamente e pode ser, por exemplo, uma determinada proporção de médias ou diferença de médias, não necessariamente padronizada).

Claramente, você já tem um tamanho de amostra, portanto não está no caso (3). Você pode considerar razoavelmente o caso (2) ou o caso (1).

Sugiro o caso (1) (que também fornece uma maneira de lidar com o caso (2)).

Para ilustrar uma abordagem do caso (1) e ver como ele se relaciona com o caso (2), vamos considerar um exemplo específico, com:

• alternativas de mudança de escala

• populações exponenciais

• tamanhos de amostra nas duas amostras de 64 e 54

Como os tamanhos das amostras são diferentes, temos que considerar o caso em que a propagação relativa em uma das amostras é menor e maior que 1 (se fossem do mesmo tamanho, considerações de simetria possibilitam considerar apenas um lado). No entanto, por estarem quase do mesmo tamanho, o efeito é muito pequeno. De qualquer forma, corrija o parâmetro para uma das amostras e varie a outra.

Então, o que se faz é:

Antecipadamente:

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

Para fazer os 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

Em R, eu fiz isso:

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 fornece a seguinte "curva" de poder

O eixo x está em uma escala de log, o eixo y é a taxa de rejeição.

É difícil dizer aqui, mas os pontos pretos são ligeiramente mais altos à esquerda do que à direita (ou seja, há um pouco mais de energia quando a amostra maior tem uma escala menor).

Usando o cdf normal inverso como uma transformação da taxa de rejeição, podemos fazer a relação entre a taxa de rejeição transformada e o log kappa (kappa está sno gráfico, mas o eixo x é escalado em log) muito quase linear (exceto próximo de 0 ) e o número de simulações foi alto o suficiente para que o ruído seja muito baixo - podemos ignorá-lo nos propósitos atuais.

Então, podemos apenas usar interpolação linear. Abaixo, são mostrados tamanhos de efeito aproximados para 50% e 80% de potência nos tamanhos de sua amostra:

Os tamanhos dos efeitos do outro lado (o grupo maior tem uma escala menor) são apenas ligeiramente alterados (podem pegar um tamanho de efeito um pouco menor), mas faz pouca diferença, então não vou entender o ponto.

Portanto, o teste detectará uma diferença substancial (de uma proporção de escalas de 1), mas não uma pequena.

Agora, para alguns comentários: não acho que os testes de hipóteses sejam particularmente relevantes para a questão de interesse subjacente ( eles são bastante semelhantes? ) E, conseqüentemente, esses cálculos de potência não nos dizem nada diretamente relevante para essa questão.

Eu acho que você aborda essa questão mais útil pré-especificando o que você acha "essencialmente o mesmo" na verdade significa operacionalmente. Isso - perseguido racionalmente a uma atividade estatística - deve levar a uma análise significativa dos dados.

Como Kolmogorov-Smirnov não é paramétrico, portanto, por definição, não pode haver análise de potência aplicável. Para ter algum tipo de estimativa, você precisa assumir um modelo de plano de fundo (e assim desviar-se do mundo não paramétrico ...) e usá-lo para calcular um dos seguintes itens: tamanho da amostra, MDE ou potência (ou seja, você conserte / escolha dois e calcule o terceiro).