Cómo simular un análisis de potencia personalizado de un modelo lm (usando R)

Siguiendo las preguntas recientes que tuvimos aquí .

Esperaba saber si alguien se había encontrado o podía compartir el código R para realizar un análisis de potencia personalizado basado en la simulación de un modelo lineal.

Más tarde, obviamente, me gustaría extenderlo a modelos más complejos, pero parece que es el lugar correcto para comenzar. Gracias.

No estoy seguro de que necesite simulación para un modelo de regresión simple. Por ejemplo, vea el documento Portable Power . Para modelos más complejos, específicamente efectos mixtos, el paquete pamm en R realiza análisis de potencia a través de simulaciones. También vea la publicación de Todd Jobe que tiene un código R para la simulación.

Aquí hay algunas fuentes de código de simulación en R. No estoy seguro si alguno aborda específicamente los modelos lineales, pero tal vez brinden un ejemplo suficiente para obtener la esencia:

• Benjamin Bolker ha escrito un gran libro datos ecológicos y modelos con R . Un borrador inicial del libro completo junto con el código Sweave está disponible en línea. El capítulo 5 aborda el análisis y la simulación de potencia.

Hay otro par de ejemplos de simulación en los siguientes sitios:

• http://www.personality-project.org/R/r.datageneration.html
• http://psy-ed.wikidot.com/simulation

Adaptado de Bolker 2009 Ecological Models and Data en R. Debe declarar la fuerza de la tendencia (es decir, la pendiente) que desea probar. Intuitivamente, una tendencia fuerte y una baja variabilidad requerirán un tamaño de muestra pequeño, una tendencia débil y una gran variabilidad requerirán un tamaño de muestra grande.

a = 2  #desired slope
b = 1  #estimated intercept
sd = 20  #estimated variability defined by standard deviation
nsim = 400  #400 simulations
pval = numeric(nsim)  #placeholder for the second for loop output
Nvec = seq(25, 100, by = 1)  #vector for the range of sample sizes to be tested
power.N = numeric(length(Nvec))   #create placeholder for first for loop output
for (j in 1:length(Nvec)) {
  N = Nvec[j]  
  x = seq(1, 20, length = Nvec[j])  #x value length needs to match sample size (Nvec) length
  for (i in 1:nsim) {   #for this value of N, create random error 400 times
    y_det = a + b * x
    y = rnorm(N, mean = y_det, sd = sd)
    m = lm(y ~ x)
    pval[i] = coef(summary(m))["x", "Pr(>|t|)"]  #all the p values for 400 sims
  }  #cycle through all N values
  power.N[j] = sum(pval < 0.05)/nsim  #the proportion of correct p-values (i.e the power)
}
power.N
plot(Nvec, power.N)  #need about 90 - 100 samples for 80% power

También puede simular cuál es la tendencia mínima que podría probar para un tamaño de muestra dado, como se muestra en el libro

bvec = seq(-2, 2, by = 0.1)
power.b = numeric(length(bvec))
for (j in 1:length(bvec)) {
  b = bvec[j]
   for (i in 1:nsim) {
     y_det = a + b * x
     y = rnorm(N, mean = y_det, sd = sd)
     m = lm(y ~ x)
     pval[i] = coef(summary(m))["x", "Pr(>|t|)"]
     }
   power.b[j] = sum(pval < 0.05)/nsim
  }
 power.b
 plot(bvec, power.b)