34 votos

Simulación de la regresión logística, el análisis de poder - los experimentos diseñados

Esta pregunta es en respuesta a una respuesta dada por @Greg Nieve en lo que respecta a una pregunta he preguntado en relación con el análisis de la potencia con la regresión logística y SAS Proc GLMPOWER.

Si estoy de diseñar un experimento y se analze los resultados en un factorial de regresión logística, ¿cómo puedo utilizar la simulación ( y aquí ) para llevar a cabo un análisis del poder?

Aquí es un simple ejemplo donde hay dos variables, la primera toma tres valores posibles {0.03, 0.06, 0.09} y el segundo es un maniquí indicador {0,1}. Para cada uno de nosotros estimamos que la tasa de respuesta para cada combinación (# de respuesta y el número de personas que comercializa a). Además, deseamos tener 3 veces la cantidad de la primera combinación de factores como los otros (que pueden ser considerados en igualdad de condiciones, debido a que esta primera combinación es nuestro probado y verdadero de la versión. Esta es una configuración como la dada en el SAS supuesto mencionado en el vinculado pregunta.

enter image description here

El modelo que se utilizará para analizar los resultados serán una regresión logística, con efectos principales y de 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 su uso con este modelo para llevar a cabo un análisis del poder?

Cuando ejecuto esto a través de SAS Proc GLMPOWER (usando STDDEV =0.05486016 que corresponde a sqrt(p(1-p)) donde p es el promedio ponderado de la muestra tasas de respuesta):

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 sólo hará uso de la clase (nominal) variables para 3, 6, 9, supra, son tratados como caracteres y que podría haber sido baja, media y alta o de otras tres cadenas. Cuando el análisis real se lleva a cabo, Var1 una numérico (y lo vamos a incluir un polinomio plazo Var1*Var1) para dar cuenta de cualquier curvatura.

La salida de SAS es

enter image description here

Así que podemos ver que necesitamos 762,112 como nuestro tamaño de la muestra (Var2 principal efecto es el más difícil de estimar) con una potencia igual a 0.80 y alfa igual a 0.05. Podríamos asignar estos a fin de que 3 veces como mucho la combinación de línea de base (es decir, 0.375 * 762112) y el resto acaba de caer igualmente en los otros 5 combinaciones.

35voto

Sean Hanley Puntos 2428

Preliminares:

  • Como se discutió en el G*Alimentación manual, hay varios tipos diferentes de poder de análisis, dependiendo de lo que desee resolver. (Es decir, $N$, el tamaño del efecto $ES$, $\alpha$, y el poder existe en relación a cada uno de los otros; especificar tres de ellos te permitirá resolver para el cuarto.)

    • en su descripción, usted quiere saber el adecuado $N$ a la captura de las tasas de respuesta especificado con $\alpha=.05$, y la potencia = 80%. Esta es un a-priori de energía.
    • podemos empezar con post-hoc de energía (determinar la potencia dada $N$, las tasas de respuesta, y alfa), ya que esto es conceptualmente más simple, y luego subir
  • En adición a @GregSnow la excelente post, otro de los grandes de la guía para la simulación de poder basadas en el análisis de los CV se puede encontrar aquí: Cálculo de la potencia estadística. Para resumir las ideas básicas:

    1. averiguar el efecto que usted quiere ser capaz de detectar
    2. generar N datos de ese mundo posible
    3. ejecute el análisis tiene la intención de llevar a cabo sobre los datos de imitación
    4. store si los resultados son "significativos" de acuerdo a su elegido alfa
    5. repetir muchas ($B$) veces y el uso de la % 'importantes' como una estimación de (post-hoc) el poder en ese $N$
    6. para determinar a priori el poder, la búsqueda de posibles $N$'s para encontrar el valor que se obtiene de su potencia deseada
  • Ya sea que usted va a encontrar significado en una determinada iteración puede ser entendido como el resultado de un ensayo de Bernoulli con una probabilidad de p $$ (donde $p$ es el poder). La proporción encontrada de más de $B$ iteraciones nos permite aproximarse a la verdad $p$. Para obtener una mejor aproximación, se puede aumentar a $B$, aunque esto también hará la simulación de tomar más tiempo.

  • En R, la principal forma de generar datos binarios con una determinada probabilidad de "éxito" es ?rbinom

    • E. g. para obtener el número de éxitos de cada 10 ensayos de Bernoulli con probabilidad p, el código sería el rbinom(n=10, size=1, prob=p), (probablemente desee asignar el resultado a una variable para el almacenamiento)
    • también puede generar este tipo de datos menos elegantemente por el uso ?runif, por ejemplo, ifelse(runif(1)<=p, 1, 0)
    • si usted cree que los resultados están mediadas por una latente de Gauss variable, se podría generar la variable latente como una función de su covariables con ?rnormy, a continuación, convertirlos en probabilidades con pnorm() y el uso de aquellos que en su rbinom() código.
  • Usted afirma que "incluyen un polinomio plazo Var1*Var1) para dar cuenta de cualquier curvatura". Hay una confusión aquí; los términos polinomiales nos puede ayudar a cuenta para la curvatura, pero este es un término de interacción, no nos ayuda en este camino. Sin embargo, sus tasas de respuesta nos obligan a incluir tanto el cuadrado de los términos y los términos de interacción en nuestro modelo. Específicamente, el modelo deberá incluir: $var1^2$, $var1*var2$, y $var1^2*var2$, más allá de los términos básicos.

  • Aunque escrito en el contexto de una pregunta, mi respuesta aquí: Diferencia entre los modelos logit y probit tiene un montón de información básica acerca de estos tipos de modelos.
  • Así como hay diferentes tipos de error de Tipo I de las tasas cuando hay varias hipótesis (por ejemplo, por el contrario la tasa de error, la tasa de error familywise, y por la familia de la tasa de error), por lo que hay diferentes tipos de poder* (por ejemplo, para una sola pre-especificado efecto, para cualquier efecto, y para todos los efectos). También puede buscar el poder para detectar una combinación específica de los efectos, o por el poder de una simultánea de prueba del modelo como un todo. Supongo que a partir de su descripción de su SAS código es que lo que busca es el último. Sin embargo, a partir de su descripción de su situación, estoy asumiendo que usted desea para detectar los efectos de la interacción en un mínimo.

  • Para una forma diferente de pensar acerca de los problemas relacionados con la alimentación, véase mi respuesta a esta pregunta: ¿Cómo reportar general de precisión en la estimación de las correlaciones dentro de un contexto de justificar el tamaño de la muestra.

Simple post-hoc de energía para la regresión logística en R:

Digamos que el postulado de las tasas de respuesta representan la verdadera situación en el mundo, y que había enviado a los 10.000 cartas. ¿Cuál es la potencia para detectar esos efectos? (Tenga en cuenta que soy famoso por escribir "cómicamente ineficiente", el código siguiente está diseñado para ser fácil de seguir, en lugar de optimizar la eficiencia; de hecho, es bastante lenta.)

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

Así que podemos ver que 10.000 cartas realmente no alcanzar el 80% de la potencia (de cualquier tipo) para detectar estas tasas de respuesta. (Yo no soy lo suficientemente seguro acerca de lo que el SAS código está haciendo para ser capaz de explicar las marcadas diferencias existentes entre estos enfoques, pero este código es conceptualmente sencillo--si lento, y yo he pasado algún tiempo, el control de ella, y creo que estos resultados son razonables.)

Basado en la simulación de un a-priori de energía para la regresión logística:

A partir de aquí, la idea es simplemente buscar una posible $N$'s hasta encontrar un valor que se obtiene el nivel deseado del tipo de alimentación que usted está interesado en. Cualquier estrategia de búsqueda que usted puede código para trabajar con este estaría bien (en teoría). Dado el $N$'s de que van a ser necesarios para capturar esos pequeños efectos, vale la pena pensar acerca de cómo hacer esto de manera más eficiente. Mi enfoque típico es simplemente la fuerza bruta, es decir, para evaluar cada $N$ que yo podría razonablemente en cuenta. (Tenga en cuenta sin embargo, que me gustaría que sólo suele considerar un rango pequeño, y estoy trabajando normalmente con muy pequeñas $N$'s--al menos en comparación con esto.)

En cambio, mi estrategia fue el de soporte posible $N$'s para obtener un sentido de lo que el rango de potencias. Por lo tanto, he elegido un $N$ 500.000 y re-corrió el código (la iniciación de la misma semilla, n.b. esto llevó a una hora y media para que se ejecute). Aquí están los resultados:

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

Aquí podemos ver que la magnitud de sus efectos varía considerablemente, y por lo tanto su capacidad para detectar ellos varía. Por ejemplo, el efecto de $var1^2$ es particularmente difícil de detectar, sólo siendo significativo el 6% del tiempo, incluso con la mitad de un millón de cartas. Por otro lado, el modelo como un todo siempre fue significativamente mejor que el modelo nulo. Las otras posibilidades son dispuestos en el medio. Aunque la mayoría de los 'datos' se tiran a la basura en cada iteración, una buena parte de la exploración es todavía posible. Por ejemplo, podríamos utilizar la significant matriz para evaluar las correlaciones entre las probabilidades de las diferentes variables que son importantes.

Debo señalar que, en conclusión, que debido a la complejidad y a la gran $N$ supuso en su situación, esto no fue tan fácil como yo había sospechado / reivindicada en mi comentario inicial. Sin embargo, sin duda se puede tener la idea de cómo esto se puede hacer en general, y de las cuestiones implicadas en el análisis del poder, de lo que he puesto aquí. HTH.

16voto

Eero Puntos 1612

@Gung la respuesta es excelente para la comprensión. Aquí es el enfoque que yo uso:

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 de las pruebas de que el efecto global de la v2, los modelos pueden ser modificados para buscar en otros tipos de pruebas. Me gusta escribir como una función de modo que cuando quiero probar algo diferente, simplemente puedo cambiar los argumentos de la función. También puede agregar una barra de progreso, o utilizar el paralelo paquete para acelerar las cosas.

Aquí acabo de hacer 100 repeticiones, por lo general comienzan en torno a ese nivel para encontrar el aproximado del tamaño de la muestra, a continuación, hasta el itterations cuando estoy en el derecho de la bola del parque (no hay necesidad de perder el tiempo en 10.000 iteraciones cuando tienes 20% de la potencia).

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by: