Análisis de potencia para regresión logística ordinal

Estoy buscando un programa (en R o SAS o independiente, si es gratuito o de bajo costo) que realizará análisis de potencia para la regresión logística ordinal.

Prefiero hacer análisis de potencia más allá de lo básico por simulación. Con los paquetes preescaneados, nunca estoy seguro de qué suposiciones se están haciendo.

Simular potencia es bastante sencillo (y asequible) con R.

1. decida cómo cree que deberían verse sus datos y cómo los analizará
2. escriba una función o un conjunto de expresiones que simularán los datos para una relación determinada y un tamaño de muestra y realice el análisis (es preferible una función en la que pueda convertir el tamaño de la muestra y los parámetros en argumentos para que sea más fácil probar diferentes valores). La función o el código deben devolver el valor p u otra estadística de prueba.
3. utilizo la replicatefunción para ejecutar el código desde arriba varias veces (generalmente comienzo aproximadamente 100 veces para tener una idea de cuánto tiempo lleva y obtener el área general correcta, luego hasta 1,000 y a veces 10,000 o 100,000 para el valores finales que usaré). La proporción de veces que rechazaste la hipótesis nula es el poder.
4. rehaga lo anterior para otro conjunto de condiciones.

Aquí hay un ejemplo simple con regresión ordinal:

library(rms)

tmpfun <- function(n, beta0, beta1, beta2) {
    x <- runif(n, 0, 10)
    eta1 <- beta0 + beta1*x
    eta2 <- eta1 + beta2
    p1 <- exp(eta1)/(1+exp(eta1))
    p2 <- exp(eta2)/(1+exp(eta2))
    tmp <- runif(n)
    y <- (tmp < p1) + (tmp < p2)
    fit <- lrm(y~x)
    fit$stats[5]
}

out <- replicate(1000, tmpfun(100, -1/2, 1/4, 1/4))
mean( out < 0.05 )

Agregaría otra cosa a la respuesta de Snow (y esto se aplica a cualquier análisis de potencia mediante simulación): preste atención a si está buscando una prueba de 1 o 2 colas. Los programas populares como G * Power tienen por defecto la prueba de 1 cola, y si está tratando de ver si sus simulaciones coinciden (siempre es una buena idea cuando está aprendiendo cómo hacer esto), querrá verificar eso primero.

Para hacer que Snow ejecute una prueba de 1 cola, agregaría un parámetro llamado "cola" a las entradas de la función, y pondría algo como esto en la función misma:

 #two-tail test
  if (tail==2) fit$stats[5]

  #one-tail test
  if (tail==1){
    if (fit$coefficients[5]>0) {
          fit$stats[5]/2
    } else 1

La versión de 1 cola básicamente verifica que el coeficiente sea positivo, y luego corta el valor p por la mitad.

Además del excelente ejemplo de Snow, creo que también puedes hacer una simulación de potencia al volver a tomar muestras de un conjunto de datos existente que tenga tu efecto. No es exactamente una rutina de arranque, ya que no estás probando con reemplazo la misma n , sino la misma idea.

Así que aquí hay un ejemplo: realicé un pequeño autoexperimento que resultó en una estimación puntual positiva, pero debido a que era pequeño, no fue estadísticamente significativo en la regresión logística ordinal. Con ese punto de presupuesto, cómo un gran n voy a necesitar? Para varias posibles n , muchas veces generé un conjunto de datos y ejecuté la regresión logística ordinal y vi cuán pequeño era el valor p :

library(boot)
library(rms)
npt <- read.csv("http://www.gwern.net/docs/nootropics/2013-gwern-noopept.csv")
newNoopeptPower <- function(dt, indices) {
    d <- dt[sample(nrow(dt), n, replace=TRUE), ] # new dataset, possibly larger than the original
    lmodel <- lrm(MP ~ Noopept + Magtein, data = d)
    return(anova(lmodel)[7])
}
alpha <- 0.05
for (n in seq(from = 300, to = 600, by = 30)) {
   bs <- boot(data=npt, statistic=newNoopeptPower, R=10000, parallel="multicore", ncpus=4)
   print(c(n, sum(bs$t<=alpha)/length(bs$t)))
}

Con la salida (para mí):

[1] 300.0000   0.1823
[1] 330.0000   0.1925
[1] 360.0000   0.2083
[1] 390.0000   0.2143
[1] 420.0000   0.2318
[1] 450.0000   0.2462
[1] 480.000   0.258
[1] 510.0000   0.2825
[1] 540.0000   0.2855
[1] 570.0000   0.3184
[1] 600.0000   0.3175

En este caso, a n = 600 la potencia fue del 32%. No muy alentador.

(Si mi enfoque de simulación es incorrecto, por favor alguien dígame. Voy a dejar algunos documentos médicos que discuten la simulación de potencia para planificar ensayos clínicos, pero no estoy del todo seguro acerca de mi implementación precisa).

En referencia a la primera simulación (sugerida por Snow; /stats//a/22410/231675 ):

Todavía no estoy seguro de cómo debería ser la simulación con más (específicamente, tres) variables independientes. Entiendo que debería 'incluirlos a todos en la pieza eta1, por ejemplo, eta1 <- beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3' '(como se mencionó anteriormente). Pero no sé cómo ajustar el resto de los parámetros en la función. ¿Alguien podría ayudarme con esto?