Simulación de análisis de potencia de regresión logística: experimentos diseñados

39

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.

ingrese la descripción de la imagen aquí

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

ingrese la descripción de la imagen aquí

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.

B_Miner
fuente
Esto es fácil de hacer en R. Primera pregunta: ¿estoy en lo cierto al decir que desea que el 75% de todos los casos sean {var1 = .03, var2 = 0} y 25% para todos los demás combos, y no 3 unidades por cada 1 unidad en cada uno de los otros combos (es decir, 37.5%)? Segunda pregunta, ¿puede especificar los efectos que le interesa detectar? Es decir, ¿cuáles serían las probabilidades de registro de 1 vs 0? ¿Cómo deberían cambiar las probabilidades de éxito del registro si var1 aumenta en .01? ¿Crees que podría haber una interacción (si es así, qué tan grande es)? (NB, estas Q pueden ser difíciles de responder, 1 estrategia es especificar la proporción de 1 que crees que podría estar en cada combo.)
gung - Restablece a Monica
Primero: El peso de 3 para el caso de referencia es que hay 3 veces más casos donde {var1 = 0.03, var2 = 0}. Por lo tanto, los resultados de SAS (que dice que necesitamos un tamaño de muestra total de 762,112 para tener un poder del 80% de rechazar el efecto principal var2 = 0, por lo que ese es el tamaño de muestra total que necesitamos) se asignaría un 37.5% a este caso de referencia.
B_Miner
2º: Bueno, todo lo que tenemos son las tasas de respuesta (que es la proporción esperada del número de éxitos sobre el número de ensayos). Entonces, si enviamos 1,000 cartas con Var1 = 0.03 y Var2 = 0, lo que podría corresponder a una oferta de tasa de interés en una oferta de correo directo de tarjeta de crédito de 0.03 (3%) y sin etiqueta en el sobre (donde Var2 = 1 significa que hay una pegatina), esperamos 1000 * 0.0025 respuestas.
B_Miner
2do cont: Esperamos una interacción, de ahí las tasas de respuesta. Tenga en cuenta que hay una tasa de respuesta diferente para Var2 = 0 dependiendo del valor de Var1. No estoy seguro de cómo traducir esto para registrar las probabilidades y luego el conjunto de datos simulados.
B_Miner
Sin embargo, una última cosa. Noto que las tasas de respuesta son lineales para var1 cuando var2 = 0 (es decir, .25%, .30%, .35%). ¿Pretendías que esto fuera un efecto lineal o curvilíneo? Debe saber que las probabilidades pueden parecer bastante lineales para pequeños subconjuntos de su rango, pero en realidad no pueden ser lineales. La regresión logística es lineal en las probabilidades de registro, no en la probabilidad (discuto cosas como esa en mi respuesta aquí ).
gung - Restablece a Monica

Respuestas:

43

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

    • 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 . norteα=.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 avanzarnorte
  • 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 esesinorte
    6. Para determinar la potencia a priori, busque posibles para encontrar el valor que le proporciona la potencia deseada. norte
  • 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 Bpagspagssipagssi

  • 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 2v a r 2vunar12vunar1vunar2vunar12vunar2

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

  • 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 Nnortenortenortenorte

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: Nnortenorte

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

vunar12significant

norte

gung - Restablece a Monica
fuente
Gung - ¡GUAU, muchas gracias por una respuesta tan detallada y reflexiva! Al escribir el mío y jugar con su código, los términos cuadráticos parecen ser el problema, ya que al menos el 80% de potencia se logra con un tamaño de muestra mucho más pequeño sin tenerlo en cuenta en el modelo.
B_Miner
1
Eso es genial, @B_Miner, ese es el tipo de cosas que quieres hacer. Además, es la razón por la que creo que el enfoque basado en la simulación es superior al software analítico que solo escupe un número (R también tiene esto, el pwrpaquete). Este enfoque le brinda la oportunidad de pensar mucho más claramente (y / o refinar su pensamiento) sobre lo que espera que suceda, cómo trataría con eso, etc. NB, sin embargo, necesita los términos cuadráticos, o algo así. análogamente, si sus tasas postuladas son correctas, b / c no son lineales, y la interacción por sí sola no le permite capturar relaciones curvilíneas.
gung - Restablece a Monica
Creo que deberías demostrar el uso de, en polylugar de mostrar a los nuevos usuarios de R, la estrategia más propensa a errores de elevar los valores sin procesar. Creo que el modelo completo debería haberse presentado como glm(responses~ poly(var1, 2) * var2, family=binomial(link="logit"). Sería menos propenso a errores estadísticos en la interpretación y mucho más compacto. Puede que no sea importante en este caso exacto cuando solo está buscando un ajuste general, pero podría confundir fácilmente a los usuarios menos sofisticados que podrían verse tentados a mirar términos individuales.
DWin
@DWin, cuando uso R para ilustrar cosas aquí en CV, lo hago de una manera muy diferente a R. La idea es ser lo más transparente posible para aquellos que no están familiarizados con R. Por ejemplo, no estoy usando las posibilidades vectorizadas, estoy usando bucles =, etc. La gente estará más familiarizada con las variables de cuadratura a partir de una regresión básica. clase, y menos conscientes de lo que poly()es, si no son usuarios de R.
gung - Restablece a Monica
17

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

Greg Snow
fuente
Gracias Greg! Me preguntaba sobre este mismo enfoque (si estoy entendiendo correctamente lo que hiciste). Para confirmar: ¿Está creando un conjunto de datos (en efecto, pero haciéndolo con pesas en lugar de fuerza bruta creando registros individuales de los valores de Var1 y Var2 y luego 1 y 0 para la respuesta) que es muy grande basado en "mydat" , ajustando una regresión logística y luego usando esos coeficientes para muestrear del modelo propuesto en la simulación? Parece que esta es una forma general de llegar a los coeficientes, entonces es como su respuesta sobre el poder de regresión ordinal al que me vinculé.
B_Miner
El modelo inicial usa pesos para obtener los coeficientes a usar, pero en la simulación está creando un marco de datos con nfilas. Puede ser más eficiente hacer pesos en la función también.
Greg Snow
¿Estoy en lo cierto al decir que está utilizando los datos inicialmente (haciendo que sea grande obtener estimaciones muy buenas) con el fin de obtener los coeficientes que se utilizan?
B_Miner
Sí, bueno, lo grande es que la proporción multiplicada por el peso da un número entero.
Greg Snow
2
@B_Miner, estoy planeando un artículo, no sé si hay suficiente para un libro completo o no. Pero gracias por el aliento.
Greg Snow