Esta pergunta é uma resposta a uma resposta dada por @Greg Snow em relação a uma pergunta que eu fiz sobre análise de potência com regressão logística e SAS Proc GLMPOWER
.
Se estou projetando um experimento e analisando os resultados em uma regressão logística fatorial, como posso usar a simulação (e aqui ) para realizar uma análise de potência?
Aqui está um exemplo simples, onde existem duas variáveis, a primeira assume três valores possíveis {0,03, 0,06, 0,09} e a segunda é um indicador fictício {0,1}. Para cada um, estimamos a taxa de resposta de cada combinação (número de respondentes / número de pessoas comercializadas). Além disso, desejamos ter três vezes mais da primeira combinação de fatores que os outros (que podem ser considerados iguais), porque essa primeira combinação é a nossa versão testada e verdadeira. Essa é uma configuração como a apresentada no curso SAS mencionado na pergunta vinculada.
O modelo que será utilizado para analisar os resultados será uma regressão logística, com principais efeitos e interação (resposta é 0 ou 1).
mod <- glm(response ~ Var1 + Var2 + I(Var1*Var2))
Como posso simular um conjunto de dados para usar com este modelo para realizar uma análise de energia?
Quando eu executo isso no SAS Proc GLMPOWER
(usando STDDEV =0.05486016
qual corresponde a sqrt(p(1-p))
onde p é a média ponderada das taxas de resposta mostradas):
data exemplar;
input Var1 $ Var2 $ response weight;
datalines;
3 0 0.0025 3
3 1 0.00395 1
6 0 0.003 1
6 1 0.0042 1
9 0 0.0035 1
9 1 0.002 1;
run;
proc glmpower data=exemplar;
weight weight;
class Var1 Var2;
model response = Var1 | Var2;
power
power=0.8
ntotal=.
stddev=0.05486016;
run;
Nota: GLMPOWER
somente usará variáveis de classe (nominais) para que 3, 6, 9 acima sejam tratadas como caracteres e poderiam ter sido baixo, médio e alto ou quaisquer outras três cadeias. Quando a análise real é realizada, Var1 será usado numérico (e incluiremos um termo polinomial Var1 * Var1) para explicar qualquer curvatura.
A saída do SAS é
Portanto, vemos que precisamos de 762.112 como o tamanho da amostra (o efeito principal da Var2 é o mais difícil de estimar) com potência igual a 0,80 e alfa igual a 0,05. Alocá-las-emos para que três vezes mais fosse a combinação da linha de base (ou seja, 0,375 * 762112) e o restante caia igualmente nas outras 5 combinações.
Respostas:
Preliminares:
Além do excelente post de @ GregSnow , outro guia realmente ótimo para análises de potência baseadas em simulação no CV pode ser encontrado aqui: Calculando o poder estatístico . Para resumir as idéias básicas:
Em R, a principal maneira de gerar dados binários com uma certa probabilidade de 'sucesso' é ? Rbinom
rbinom(n=10, size=1, prob=p)
(você provavelmente desejará atribuir o resultado a uma variável para armazenamento)ifelse(runif(1)<=p, 1, 0)
pnorm()
e usá-las no seurbinom()
código.Assim como existem diferentes tipos de taxas de erro tipo I quando há várias hipóteses (por exemplo, taxa de erro por contraste , a taxa de erro familywise , e taxa de erro por família ), de modo que existem diferentes tipos de poder * (por exemplo, para um único efeito pré-especificado , para qualquer efeito e para todos os efeitos ). Você também pode procurar o poder de detectar uma combinação específica de efeitos ou o poder de um teste simultâneo do modelo como um todo. Meu palpite, pela descrição do seu código SAS, é que ele está procurando pelo último. No entanto, a partir da descrição da sua situação, suponho que você queira detectar os efeitos da interação no mínimo.
Para uma maneira diferente de pensar sobre questões relacionadas ao poder, veja minha resposta aqui: Como relatar precisão geral na estimativa de correlações em um contexto de justificação do tamanho da amostra.
Poder post-hoc simples para regressão logística em R:
Digamos que suas taxas de resposta representam a verdadeira situação no mundo e que você enviou 10.000 cartas. Qual é o poder de detectar esses efeitos? (Observe que eu sou famoso por escrever código "comicamente ineficiente", o seguinte deve ser fácil de seguir, em vez de otimizado para eficiência; na verdade, é bem lento.)
Portanto, vemos que 10.000 letras não atingem 80% de energia (de qualquer tipo) para detectar essas taxas de resposta. (Não tenho certeza suficiente sobre o que o código SAS está fazendo para poder explicar a grande discrepância entre essas abordagens, mas esse código é conceitualmente direto - se lento - e passei algum tempo verificando-o e acho que esses os resultados são razoáveis.)
Poder a priori baseado em simulação para regressão logística:
significant
fonte
pwr
pacote). Essa abordagem oferece a oportunidade de pensar com muito mais clareza (e / ou refinar seu pensamento) sobre o que você espera que aconteça, como você lida com isso, etc. NB, no entanto, você precisa dos termos quadráticos, ou algo assim de forma análoga, se suas taxas postadas estiverem corretas, b / c elas não são lineares e a interação por si só não permite capturar relacionamentos curvilíneos.poly
vez de mostrar aos novos usuários de R, a estratégia mais propensa a erros de quadratura de valores brutos. Eu acho que o modelo completo deveria ter sido colocado comoglm(responses~ poly(var1, 2) * var2, family=binomial(link="logit")
. Seria menos propenso a erros estatísticos na interpretação e muito mais compacto. Pode não ser importante neste caso exato quando você está apenas olhando para um ajuste geral, mas pode facilmente enganar usuários menos sofisticados que podem ficar tentados a olhar para termos individuais.=
, etc. As pessoas estarão mais familiarizadas com variáveis quadráticas a partir de uma regressão básica classe e menos ciente do quepoly()
é, se não forem usuários R.A resposta de @ Gung é ótima para entender. Aqui está a abordagem que eu usaria:
Esta função testa o efeito geral da v2, os modelos podem ser alterados para observar outros tipos de testes. Eu gosto de escrevê-lo como uma função para que, quando eu quiser testar algo diferente, eu possa mudar os argumentos da função. Você também pode adicionar uma barra de progresso ou usar o pacote paralelo para acelerar as coisas.
Aqui, eu fiz apenas 100 repetições, normalmente começo nesse nível para encontrar o tamanho aproximado da amostra e, em seguida, subo as iterações quando estou no parque de bola certo (não há necessidade de perder tempo em 10.000 iterações quando você tem 20% de energia).
fonte
n
linhas. Pode ser mais eficiente fazer pesos na função também.