Esta pregunta es en respuesta a una respuesta dada por @Greg Snow con respecto a una pregunta que hice sobre el análisis de potencia con regresión logística y SAS Proc GLMPOWER.

Si estoy diseñando un experimento y analizaré los resultados en una regresión logística factorial, ¿cómo puedo usar la simulación (y aquí ) para realizar un análisis de potencia?

Aquí hay un ejemplo simple donde hay dos variables, la primera toma tres valores posibles {0.03, 0.06, 0.09} y la segunda es un indicador ficticio {0,1}. Para cada uno, estimamos la tasa de respuesta para cada combinación (# de respondedores / número de personas comercializadas). Además, deseamos tener 3 veces más de la primera combinación de factores que los otros (que puede considerarse igual) porque esta primera combinación es nuestra versión probada y verdadera. Esta es una configuración como la dada en el curso SAS mencionado en la pregunta vinculada.

El modelo que se utilizará para analizar los resultados será una regresión logística, con efectos principales e interacción (la respuesta es 0 o 1).

mod <- glm(response ~ Var1 + Var2 + I(Var1*Var2))

¿Cómo puedo simular un conjunto de datos para usar con este modelo para realizar un análisis de potencia?

Cuando ejecuto esto a través de SAS Proc GLMPOWER(usando STDDEV =0.05486016 cuál corresponde a sqrt(p(1-p))donde p es el promedio ponderado de las tasas de respuesta 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: GLMPOWERsolo se usarán variables de clase (nominales), por lo que 3, 6, 9 anteriores se tratan como caracteres y podrían haber sido bajos, medios y altos o cualquier otra cadena. Cuando se realiza el análisis real, Var1 se utilizará de forma numérica (e incluiremos un término polinómico Var1 * Var1) para dar cuenta de cualquier curvatura.

La salida de SAS es

Entonces vemos que necesitamos 762,112 como nuestro tamaño de muestra (el efecto principal de Var2 es el más difícil de estimar) con una potencia igual a 0,80 y alfa igual a 0,05. Los asignaríamos de modo que 3 veces más fueran la combinación de línea de base (es decir, 0.375 * 762112) y el resto caiga igualmente en las otras 5 combinaciones.

Preliminares:

• Como se discutió en el manual de G * Power , hay varios tipos diferentes de análisis de potencia, dependiendo de lo que desee resolver. (Es decir, , el tamaño del efecto , y la potencia existen entre sí; especificar tres de ellos te permitirá resolver el cuarto). $N$$ES$$\alpha$

  • en su descripción, desea conocer el apropiado para capturar las tasas de respuesta que especificó con y power = 80%. Este es un poder a priori . $N$$\alpha=.05$
  • podemos comenzar con la potencia post-hoc (determinar la potencia dada , las tasas de respuesta y alfa) ya que esto es conceptualmente más simple, y luego avanzar$N$

• Además de la excelente publicación de @ GregSnow, aquí se puede encontrar otra guía realmente excelente para análisis de potencia basados ​​en simulación en CV: Cálculo de potencia estadística . Para resumir las ideas básicas:

  1. averiguar el efecto que desea poder detectar
  2. generar N datos de ese mundo posible
  3. ejecuta el análisis que pretendes realizar sobre esos datos falsos
  4. almacenar si los resultados son 'significativos' de acuerdo con el alfa elegido
  5. repetir muchas ( ) veces y usar el% 'significativo' como una estimación de la potencia (post-hoc) en ese$B$$N$
  6. Para determinar la potencia a priori, busque posibles para encontrar el valor que le proporciona la potencia deseada. $N$

• Si usted encontrará importancia en una iteración particular puede entenderse como el resultado de un ensayo de Bernoulli con probabilidad (donde es el poder). La proporción encontrada sobre iteraciones nos permite aproximar la verdadera . Para obtener una mejor aproximación, podemos aumentar , aunque esto también hará que la simulación tarde más. p B p $p$$p$$B$$p$$B$

• En R, la forma principal de generar datos binarios con una probabilidad dada de 'éxito' es ? Rbinom

  • Por ejemplo, para obtener el número de éxitos de 10 ensayos de Bernoulli con probabilidad p, el código sería rbinom(n=10, size=1, prob=p)(probablemente querrá asignar el resultado a una variable para el almacenamiento)
  • También puede generar dichos datos de manera menos elegante utilizando ? runif , por ejemplo,ifelse(runif(1)<=p, 1, 0)
  • si cree que los resultados están mediados por una variable Gaussiana latente, podría generar la variable latente en función de sus covariables con ? rnorm , y luego convertirlas en probabilidades pnorm()y utilizarlas en su rbinom()código.

• Usted declara que "incluirá un término polinómico Var1 * Var1) para dar cuenta de cualquier curvatura". Hay una confusión aquí; Los términos polinómicos pueden ayudarnos a explicar la curvatura, pero este es un término de interacción, no nos ayudará de esta manera. No obstante, sus tasas de respuesta requieren que incluyamos términos cuadrados y términos de interacción en nuestro modelo. Específicamente, su modelo deberá incluir: , y , más allá de los términos básicos. v a r 1 ∗ v a r 2 v a r 1 2 ∗ v a r $var1^2$$var1*var2$$var1^2*var2$

• Aunque escrito en el contexto de una pregunta diferente, mi respuesta aquí: la diferencia entre los modelos logit y probit tiene mucha información básica sobre este tipo de modelos.

• Al igual que hay diferentes tipos de tasas de error tipo I cuando hay varias hipótesis (por ejemplo, tasa de error por contraste , la tasa de errores familywise , y tasa de error por familia ), por lo que hay diferentes tipos de poder * (por ejemplo, para una solo efecto preespecificado , para cualquier efecto y para todos los efectos ). También puede buscar el poder de detectar una combinación específica de efectos, o el poder de una prueba simultánea del modelo en su conjunto. Supongo que por su descripción de su código SAS es que está buscando el último. Sin embargo, a partir de su descripción de su situación, supongo que desea detectar los efectos de interacción como mínimo.

  • * referencia: Maxwell, SE (2004). La persistencia de estudios de baja potencia en investigación psicológica: causas, consecuencias y remedios. Métodos psicológicos , 9 , 2 , pp. 147-163.
  • sus efectos son bastante pequeños (no debe confundirse con las bajas tasas de respuesta), por lo que nos será difícil lograr un buen poder.
  • Tenga en cuenta que, aunque todo esto suena bastante similar, en gran medida no son lo mismo (por ejemplo, es muy posible obtener un modelo significativo sin efectos significativos - discutido aquí: ¿Cómo puede una regresión ser significativa, pero todos los predictores son no- significativo?, o efectos significativos pero donde el modelo no es significativo - discutido aquí: Importancia de los coeficientes en regresión lineal: prueba t significativa versus estadística F no significativa ), que se ilustrará a continuación.

• Para ver una forma diferente de pensar sobre cuestiones relacionadas con el poder, vea mi respuesta aquí: Cómo informar la precisión general en la estimación de correlaciones dentro de un contexto de justificación del tamaño de la muestra.

Potencia post-hoc simple para regresión logística en R:

Digamos que sus tasas de respuesta postuladas representan la verdadera situación en el mundo, y que ha enviado 10,000 cartas. ¿Cuál es el poder de detectar esos efectos? (Tenga en cuenta que soy famoso por escribir código "cómicamente ineficiente", lo siguiente pretende ser fácil de seguir en lugar de estar optimizado para la eficiencia; de hecho, es bastante 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

Entonces, vemos que 10,000 letras realmente no alcanzan el 80% de potencia (de ningún tipo) para detectar estas tasas de respuesta. (No estoy lo suficientemente seguro sobre lo que está haciendo el código SAS para poder explicar la clara discrepancia entre estos enfoques, pero este código es conceptualmente directo, si es lento, y he pasado algún tiempo revisándolo, y creo que estos los resultados son razonables)

Potencia a priori basada en simulación para regresión logística:

A partir de aquí, la idea es simplemente buscar posibles 's hasta que encontremos un valor que produzca el nivel deseado del tipo de poder que le interesa. Cualquier estrategia de búsqueda que pueda codificar para trabajar con esto estaría bien (en teoría). Dadas las que se necesitarán para capturar efectos tan pequeños, vale la pena pensar en cómo hacer esto de manera más eficiente. Mi enfoque típico es simplemente la fuerza bruta, es decir, evaluar cada que podría considerar razonablemente. (Sin embargo, tenga en cuenta que normalmente solo consideraría un rango pequeño, y normalmente estoy trabajando con muy pequeñas , al menos en comparación con esto). N N $N$$N$$N$$N$

En cambio, mi estrategia aquí fue poner entre paréntesis posibles para tener una idea de cuál sería el rango de poderes. Por lo tanto, elegí un de 500,000 y volví a ejecutar el código (iniciando la misma semilla, esto tardó una hora y media en ejecutarse). Aquí están los resultados: $N$$N$

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

$var1^2$significant

$N$

La respuesta de @ Gung es genial para entender. Aquí está el enfoque que usaría:

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 función prueba el efecto general de v2, los modelos se pueden cambiar para ver otros tipos de pruebas. Me gusta escribirlo como una función para que cuando quiera probar algo diferente pueda cambiar los argumentos de la función. También puede agregar una barra de progreso o usar el paquete paralelo para acelerar las cosas.

Aquí acabo de hacer 100 repeticiones, generalmente comienzo alrededor de ese nivel para encontrar el tamaño aproximado de la muestra, y luego las iteraciones cuando estoy en el parque de pelota correcto (no es necesario perder el tiempo en 10,000 iteraciones cuando tienes un 20% de potencia).