Cálculo do tamanho da amostra para regressão logística univariada

11

Como se calcula o tamanho da amostra necessário para um estudo no qual uma coorte de sujeitos terá uma única variável contínua medida no momento da cirurgia e, dois anos depois, eles serão classificados como resultado funcional ou resultado prejudicado.

Gostaríamos de ver se essa medida poderia ter previsto o resultado ruim. Em algum momento, podemos derivar um ponto de corte na variável contínua acima da qual tentaríamos intervir para diminuir a probabilidade do resultado prejudicado.

Alguma ideia? Qualquer implementação de R.

Farrel
fonte
Você espera algumas desistências durante o acompanhamento? Existem outras covariáveis ​​a serem incluídas no seu modelo?
chl
Deixe-me sugar uma taxa de evasão do meu polegar - 20%. De fato, coletaremos muitas variáveis, por exemplo, idade, pontuação no trauma, mas eu queria manter as coisas o mais simples possível para o cálculo da potência. Muitas vezes achei útil discutir um modelo primário e depois modelos secundários carregados com mais requinte e nuances.
Farrel
Ok, mas geralmente a porcentagem de evasão esperada, o número de covariáveis ​​e se as covariáveis ​​são medidas com erros (veja por exemplo, j.mp/9fJkhb ) inserem a fórmula (em todos os casos, isso aumentará o tamanho da amostra).
chl

Respostas:

7

Os cálculos do tamanho da amostra para regressão logística são complexos. Não vou tentar resumir aqui. Soluções razoavelmente acessíveis para esse problema são encontradas em:

Hsieh FY. Tabelas de tamanho de amostra para regressão logística. Estatística em Medicina. Julho de 1989; 8 (7): 795-802.

Hsieh FY, et al. Um método simples de cálculo do tamanho da amostra para regressão linear e logística. Estatística em Medicina. 30 de julho de 1998; 17 (14): 1623-34.

Uma discussão acessível dos problemas com exemplos de cálculos pode ser encontrada no último capítulo (Seção 8.5, pp. 339-347) da Regressão logística aplicada de Hosmer & Lemeshow .

Thylacoleo
fonte
7

Normalmente, acho mais fácil e rápido executar uma simulação. Os artigos demoram muito tempo para ler, entender e finalmente chegar à conclusão de que não se aplicam no caso especial em que estamos interessados.

Portanto, eu apenas escolhia uma série de assuntos, simulava a covariável em que você está interessado (distribuído como você acredita que seria), simula resultados bons / ruins com base na forma funcional que você posiciona (efeitos de limiar da covariável? Não linearidade?) com o tamanho mínimo (clinicamente) de efeito significativo que você gostaria de detectar, execute o resultado através de sua análise e veja se o efeito é encontrado em seu alfa. Execute novamente 10.000 vezes e verifique se você encontrou o efeito em 80% das simulações (ou qualquer outro poder que você precise). Ajuste o número de assuntos, repita até ter um poder que lhe agrade.

Isso tem a vantagem de ser muito geral, portanto você não se limita a uma forma funcional específica ou a um número ou distribuição específica de covariáveis. Você pode incluir desistências, veja o comentário de chl acima, aleatoriamente ou influenciado por covariáveis ​​ou resultados. Você basicamente codifica a análise que fará na amostra final antecipadamente, o que às vezes ajuda a focar meu pensamento no desenho do estudo. E isso é feito facilmente em R (vetorizar!).

Stephan Kolassa
fonte
Você tem um caso resolvido em R?
Farrel
1
@ Farrel - aqui está um script muito curto, que assume covariáveis ​​distribuídas [0,1] uniformemente, um OR de 2 entre o primeiro e o terceiro quartil do covariado e o ruído normal padrão, levando à potência 0,34 para n = 100. Eu brincaria com isso para ver como tudo é sensível às minhas suposições: corre <- 1000; nn <- 100; set.seed (2010); detecções <- replicar (n = execuções, expr = {covariável <- runif (nn); resultado <- runif (nn) <1 / (1 + exp (-2 * log (2) * covariável + rnorm (nn)) ); resumo (glm (resultado ~ covariável, família = "binomial")) $ coeficientes ["covariável", "Pr (> | z |)"] <.05}) cat ("Potência:", soma (detecções) / execuções, "\ n")
Stephan Kolassa
1
Você pode anexar seu código como um pastie ( pastebin.com ) ou um Gist ( gist.github.com ) se achar que é mais conveniente e vincular a ele no seu comentário.
chl
@chl: +1, muito obrigado! Aqui está a essência: gist.github.com/607968
Stephan Kolassa
Ótimo código, mas há um problema. Eu não sou tão inteligente quanto você. Eu preciso disso dividido em etapas. Eu acho que é executado o número de simulações? O que é nn? É o número de sujeitos no estudo? Então, vejo que você criou uma distribuição de covariáveis ​​e as fez determinar um sim ou um não, dependendo de um limite.
Farrel
4

Seguindo o post de Stephan Kolassa (não posso adicionar isso como comentário), tenho um código alternativo para uma simulação. Isso usa a mesma estrutura básica, mas é expandida um pouco mais, então talvez seja um pouco mais fácil de ler. Também é baseado no código de Kleinman e Horton para simular a regressão logística.

nn é o número na amostra. A covariável deve ser continuamente distribuída normalmente e padronizada para significar 0 e sd 1. Usamos rnorm (nn) para gerar isso. Selecionamos um odds ratio e o armazenamos em odds.ratio. Também escolhemos um número para a interceptação. A escolha desse número controla qual proporção da amostra experimenta o "evento" (por exemplo, 0,1, 0,4, 0,5). Você tem que brincar com esse número até obter a proporção certa. O código a seguir fornece uma proporção de 0,1 com um tamanho de amostra de 950 e um OR de 1,5:

nn <- 950
runs <- 10000
intercept <- log(9)
odds.ratio <- 1.5
beta <- log(odds.ratio)
proportion  <-  replicate(
              n = runs,
              expr = {
                  xtest <- rnorm(nn)
                  linpred <- intercept + (xtest * beta)
                  prob <- exp(linpred)/(1 + exp(linpred))
                  runis <- runif(length(xtest),0,1)
                  ytest <- ifelse(runis < prob,1,0)
                  prop <- length(which(ytest <= 0.5))/length(ytest)
                  }
            )
summary(proportion)

O resumo (proporção) confirma que a proporção é ~ 0,1

Em seguida, usando as mesmas variáveis, a potência é calculada em 10.000 execuções:

result <-  replicate(
              n = runs,
              expr = {
                  xtest <- rnorm(nn)
                  linpred <- intercept + (xtest * beta)
                  prob <- exp(linpred)/(1 + exp(linpred))
                  runis <- runif(length(xtest),0,1)
                  ytest <- ifelse(runis < prob,1,0)
                  summary(model <- glm(ytest ~ xtest,  family = "binomial"))$coefficients[2,4] < .05
                  }
            )
print(sum(result)/runs)

Eu acho que esse código está correto - eu comparei com os exemplos dados em Hsieh, 1998 (tabela 2), e parece concordar com os três exemplos dados lá. Eu também o testei contra o exemplo nas páginas 342 - 343 de Hosmer e Lemeshow, onde encontrou uma potência de 0,75 (em comparação com 0,8 em Hosmer e Lemeshow). Portanto, em algumas circunstâncias, essa abordagem subestima o poder. No entanto, quando executei o mesmo exemplo nesta calculadora on-line , descobri que ela concorda comigo e não com o resultado em Hosmer e Lemeshow.

Se alguém puder nos dizer por que esse é o caso, eu estaria interessado em saber.

Andrew
fonte
Eu tenho duas perguntas se você não se importa. 1) A função de proporção é simplesmente para obter a interceptação correta? 2) qual é a lógica por trás do uso do ytest (comparando prob a um sorteio aleatório)?
21811 B_Miner
@B_Miner 1) Ao contrário - para obter a proporção correta, você precisa definir a interceptação corretamente - então ajuste a interceptação até obter a proporção que você espera. 2) A lógica do ytest é que precisamos obter um resultado dicotômico 0 ou 1. Então, comparamos cada amostra da distribuição uniforme com a probabilidade (prob) para obter nosso resultado dicotômico. Os 'runis' não precisam ser extraídos da distribuição uniforme aleatória - um binômio ou outra distribuição pode fazer mais sentido para seus dados. Espero que isso ajude (desculpe pelo atraso na resposta).
23411 Andrew Andrew
3

Uma pergunta simples sobre o tamanho da amostra é: qual o tamanho de uma amostra necessária para obter um intervalo de confiança de 95% não superior a 2d para a média [desconhecida] da distribuição de dados. outra variante é: quão grande é uma amostra para ter potência 0,9 em ao testar H . você não parece especificar nenhum critério para escolher um tamanho de amostra.0 : θ = 0θ=10:θ=0

na verdade, parece que seu estudo será conduzido de forma seqüencial. nesse caso, pode pagar para tornar isso uma parte explícita do experimento. a amostragem seqüencial geralmente pode ser mais eficiente do que um experimento fixo de tamanho de amostra [menos observações necessárias, em média].

farrel: estou adicionando isso em resposta ao seu comentário.

para obter um tamanho de amostra, geralmente se especifica algum tipo de critério de precisão para que uma estimativa [como comprimento de um IC] OU potência em uma alternativa especificada de um teste seja executada nos dados. você parece ter mencionado esses dois critérios. não há nada de errado com isso, em princípio: basta fazer dois cálculos de tamanho de amostra - um para alcançar a precisão de estimativa desejada - e outro para obter a potência desejada na alternativa declarada. então, o maior dos dois tamanhos de amostra é o necessário. [além de dizer 80% da potência - você não mencionou o teste que planeja executar - ou a alternativa na qual deseja a potência de 80%.]

quanto ao uso da análise sequencial: se os participantes estão incluídos no estudo ao mesmo tempo, faz sentido um tamanho fixo da amostra. mas se os assuntos forem poucos e distantes, pode levar um ano ou dois [ou mais] para obter o número necessário. assim, o julgamento poderia durar três ou quatro anos [ou mais]. nesse caso, um esquema seqüencial oferece a possibilidade de parar mais cedo do que isso - se o (s) efeito (s) que você está procurando se tornar estatisticamente significativo no início do teste.

ronaf
fonte
Os critérios terão 10% de diferença na probabilidade de resultados bons ou ruins. Ou, digamos, como será regressão logística, razão de chances = 2. alfa = 0,05, potência = 80%, ainda não sei qual é a variação combinada da variável contínua, mas vamos assumir que o desvio padrão é de 7 mmHg. A análise seqüencial seria boa, mas o resultado final é de dois anos após a medição.
Farrel