Simulação de análise de poder de regressão logística - experimentos projetados

39

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.

insira a descrição da imagem aqui

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: GLMPOWERsomente 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 é

insira a descrição da imagem aqui

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.

B_Miner
fonte
Isso é fácil de fazer na R. 1ª pergunta: estou certo de que você deseja que 75% de todos os casos sejam {var1 = 0,03, var2 = 0} e 25% para todos os outros combos, e não 3 unidades lá para cada 1 unidade em cada um dos outros combos (ou seja, 37,5%)? 2ª pergunta, você pode especificar os efeitos que está interessado em detectar? Ou seja, quais seriam as chances de log de 1 vs 0? Como as chances de sucesso do log devem mudar se var1 aumenta em 0,01? Você acha que pode haver uma interação (se sim, qual o tamanho)? (NB, esses Q's podem ser difíceis de responder, 1 estratégia é especificar a proporção de 1s que você acha que pode estar em cada combinação.)
gung - Reinstate Monica
1º: o peso de 3 para o caso da linha de base é que há 3 vezes mais casos em que {var1 = 0,03, var2 = 0}. Portanto, os resultados do SAS (que diz que precisamos de 762.112 tamanho total da amostra para ter 80% de poder de rejeitar o efeito principal var2 = 0, de modo que é o tamanho total da amostra que precisamos) seriam alocados 37,5% para esse caso de linha de base.
B_Miner 9/09/12
2º: Bem, tudo o que temos são as taxas de resposta (que é a taxa esperada entre o número de sucessos e o número de tentativas). Portanto, se enviarmos 1.000 cartas com Var1 = 0,03 e Var2 = 0, o que poderia corresponder a uma oferta de taxa de juros em uma oferta de mala direta de cartão de crédito de 0,03 (3%) e sem adesivo no envelope (onde Var2 = 1 significa que há um adesivo), esperamos 1000 * 0,0025 respostas.
B_Miner 9/09/12
2º cont: Esperamos uma interação - daí as taxas de resposta. Observe que há uma taxa de resposta diferente para Var2 = 0, dependendo do valor de Var1. Não sei como traduzi-las para registrar probabilidades e, em seguida, para um conjunto de dados simulados.
B_Miner 9/09/12
Uma última coisa, no entanto. Percebo que as taxas de resposta são lineares para var1 quando var2 = 0 (ou seja, 0,25%, 0,30%, 0,35%). Você pretendia que isso fosse um efeito linear ou curvilíneo? Você deve saber que as probabilidades podem parecer razoavelmente lineares para pequenos subconjuntos de seu intervalo, mas na verdade não podem ser lineares. A regressão logística é linear nas probabilidades de log, não na probabilidade (discuto coisas assim na minha resposta aqui ).
gung - Restabelece Monica

Respostas:

43

Preliminares:

  • NESα

    • Nα=.05
    • N
  • 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:

    1. descobrir o efeito que você deseja detectar
    2. gerar N dados desse mundo possível
    3. execute a análise que você pretende realizar sobre esses dados falsos
    4. armazene se os resultados são 'significativos' de acordo com o alfa escolhido
    5. BN
    6. N
  • ppBpB

  • Em R, a principal maneira de gerar dados binários com uma certa probabilidade de 'sucesso' é ? Rbinom

    • Por exemplo, para obter o número de sucessos em 10 tentativas de Bernoulli com probabilidade p, o código seria rbinom(n=10, size=1, prob=p)(você provavelmente desejará atribuir o resultado a uma variável para armazenamento)
    • você também pode gerar esses dados com menos elegância usando ? runif , por exemplo,ifelse(runif(1)<=p, 1, 0)
    • se você acredita que os resultados são mediados por uma variável gaussiana latente, você pode gerar a variável latente em função de suas covariáveis ​​com ? rnorm e depois convertê-las em probabilidades pnorm()e usá-las no seu rbinom()código.
  • var12vumar1vumar2vumar12vumar2

  • Embora tenha sido escrita no contexto de uma pergunta diferente, minha resposta aqui: A diferença entre os modelos logit e probit possui muitas informações básicas sobre esses tipos de modelos.
  • 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.)

set.seed(1)

repetitions = 1000
N = 10000
n = N/8
var1  = c(   .03,    .03,    .03,    .03,    .06,    .06,    .09,   .09)
var2  = c(     0,      0,      0,      1,      0,      1,      0,     1)
rates = c(0.0025, 0.0025, 0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002)

var1    = rep(var1, times=n)
var2    = rep(var2, times=n)
var12   = var1**2
var1x2  = var1 *var2
var12x2 = var12*var2

significant = matrix(nrow=repetitions, ncol=7)

startT = proc.time()[3]
for(i in 1:repetitions){
  responses          = rbinom(n=N, size=1, prob=rates)
  model              = glm(responses~var1+var2+var12+var1x2+var12x2, 
                           family=binomial(link="logit"))
  significant[i,1:5] = (summary(model)$coefficients[2:6,4]<.05)
  significant[i,6]   = sum(significant[i,1:5])
  modelDev           = model$null.deviance-model$deviance
  significant[i,7]   = (1-pchisq(modelDev, 5))<.05
}
endT = proc.time()[3]
endT-startT

sum(significant[,1])/repetitions      # pre-specified effect power for var1
[1] 0.042
sum(significant[,2])/repetitions      # pre-specified effect power for var2
[1] 0.017
sum(significant[,3])/repetitions      # pre-specified effect power for var12
[1] 0.035
sum(significant[,4])/repetitions      # pre-specified effect power for var1X2
[1] 0.019
sum(significant[,5])/repetitions      # pre-specified effect power for var12X2
[1] 0.022
sum(significant[,7])/repetitions      # power for likelihood ratio test of model
[1] 0.168
sum(significant[,6]==5)/repetitions   # all effects power
[1] 0.001
sum(significant[,6]>0)/repetitions    # any effect power
[1] 0.065
sum(significant[,4]&significant[,5])/repetitions   # power for interaction terms
[1] 0.017

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:

NNNN

NN

sum(significant[,1])/repetitions      # pre-specified effect power for var1
[1] 0.115
sum(significant[,2])/repetitions      # pre-specified effect power for var2
[1] 0.091
sum(significant[,3])/repetitions      # pre-specified effect power for var12
[1] 0.059
sum(significant[,4])/repetitions      # pre-specified effect power for var1X2
[1] 0.606
sum(significant[,5])/repetitions      # pre-specified effect power for var12X2
[1] 0.913
sum(significant[,7])/repetitions      # power for likelihood ratio test of model
[1] 1
sum(significant[,6]==5)/repetitions   # all effects power
[1] 0.005
sum(significant[,6]>0)/repetitions    # any effect power
[1] 0.96
sum(significant[,4]&significant[,5])/repetitions   # power for interaction terms
[1] 0.606

vumar12significant

N

- Reinstate Monica
fonte
Gung - WOW, muito obrigado por uma resposta tão detalhada e atenciosa! Ao escrever o meu próprio e brincar com o seu código, os termos quadráticos parecem ser o problema - pois pelo menos 80% da energia é alcançada com um tamanho de amostra muito menor sem considerá-lo no modelo.
B_Miner 10/09/12
1
Isso é ótimo, @B_Miner, esse é o tipo de coisa que você deseja fazer. Além disso, é a razão pela qual acho que a abordagem baseada em simulação é superior ao software analítico que apenas expõe um número (R também possui isso, o pwrpacote). 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.
gung - Restabelece Monica
Eu acho que você deveria demonstrar o uso de, em polyvez 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 como glm(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.
DWin
@ DWin, quando eu uso R para ilustrar as coisas aqui no CV, faço isso de uma maneira muito não-R. A idéia é ser o mais transparente possível para aqueles que não estão familiarizados com R. Por exemplo, eu não estou usando as possibilidades vetorizadas, estou usando loops =, 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 que poly()é, se não forem usuários R.
gung - Restabelece Monica
17

A resposta de @ Gung é ótima para entender. Aqui está a abordagem que eu usaria:

mydat <- data.frame( v1 = rep( c(3,6,9), each=2 ),
    v2 = rep( 0:1, 3 ), 
    resp=c(0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002) )

fit0 <- glm( resp ~ poly(v1, 2, raw=TRUE)*v2, data=mydat,
    weight=rep(100000,6), family=binomial)
b0 <- coef(fit0)


simfunc <- function( beta=b0, n=10000 ) {
    w <- sample(1:6, n, replace=TRUE, prob=c(3, rep(1,5)))
    mydat2 <- mydat[w, 1:2]
    eta <- with(mydat2,  cbind( 1, v1, 
                v1^2, v2,
                v1*v2,
                v1^2*v2 ) %*% beta )
    p <- exp(eta)/(1+exp(eta))
    mydat2$resp <- rbinom(n, 1, p)

    fit1 <- glm( resp ~ poly(v1, 2)*v2, data=mydat2,
        family=binomial)
    fit2 <- update(fit1, .~ poly(v1,2) )
    anova(fit1,fit2, test='Chisq')[2,5]
}

out <- replicate(100, simfunc(b0, 10000))
mean( out <= 0.05 )
hist(out)
abline(v=0.05, col='lightgrey')

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).

Greg Snow
fonte
Obrigado Greg! Eu estava pensando sobre essa mesma abordagem (se estou entendendo corretamente o que você fez). Para confirmar: Você está criando um conjunto de dados (com efeito, mas fazendo isso com pesos, em vez de força bruta, criando registros individuais dos valores de Var1 e Var2 e, em seguida, 1 e 0 para a resposta) muito grande com base em "mydat" , ajustando uma regressão logística e depois usando esses coeficientes para amostrar a partir do modelo proposto na simulação? Parece que esta é uma maneira geral de apresentar os coeficientes - então é exatamente como a sua resposta sobre o poder de regressão ordinal ao qual vinculei.
B_Miner 10/09/12
O modelo inicial usa pesos para obter os coeficientes a serem usados, mas na simulação está criando um quadro de dados com nlinhas. Pode ser mais eficiente fazer pesos na função também.
Greg Snow
Estou certo de que você está usando os dados inicialmente (aumentando para obter estimativas muito boas) com o objetivo de obter os coeficientes usados?
B_Miner 10/09/12
Sim, bem, o grande é que a proporção vezes o peso dê um número inteiro.
Greg Snow
2
@ B_Miner, estou planejando um artigo, não sei se há o suficiente para um livro completo ou não. Mas obrigado pelo incentivo.
Greg Neve