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:
- averiguar el efecto que usted quiere ser capaz de detectar
- generar N datos de ese mundo posible
- ejecute el análisis tiene la intención de llevar a cabo sobre los datos de imitación
- store si los resultados son "significativos" de acuerdo a su elegido alfa
- repetir muchas ($B$) veces y el uso de la % 'importantes' como una estimación de (post-hoc) el poder en ese $N$
- 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.
- *referencia: Maxwell, S.E. (2004). La persistencia de estudios limitados en la investigación psicológica: causas, consecuencias y remedios. Métodos psicológicos, 9, 2, pp 147-163.
- sus efectos son muy pequeñas (que no debe confundirse con la baja tasa de respuesta), por lo que será difícil lograr una buena alimentación.
- Tenga en cuenta que, a pesar de todos estos sonido bastante similar, son muy no es el mismo (por ejemplo, es muy posible obtener un modelo significativo con efectos significativos--discutido aquí: ¿Cómo puede una regresión de ser significativo, sin embargo, de todos los predictores no es significativa?, o efectos significativos pero donde el modelo no es significativo--discutido aquí: la Significación de los coeficientes de la regresión lineal: un t-test vs no significativo de la estadística F), que se ilustra a continuación.
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.
