El poder de análisis de regresión logística ordinal

Question

Estoy buscando un programa (en el R o SAS o independiente, si es gratuito o de bajo costo) que va a hacer el análisis del poder para ordinal de regresión logística.

Answer 1

Yo prefiero poder hacer análisis más allá de los conceptos básicos de la simulación. Con prefabricados de paquetes, nunca estoy muy seguro de qué suposiciones son.

La simulación por el poder es muy sencillo (y económico) mediante R.

1. decida lo que usted piensa que sus datos debería ser y cómo va a analizar
2. escribir una función o un conjunto de expresiones que se va a simular los datos para una determinada relación y tamaño de la muestra y realizar el análisis (una función que es preferible, en la que puede hacer que el tamaño de la muestra y los parámetros en los argumentos para hacer más fácil probar diferentes valores). La función o el código debe devolver el valor de la p o la otra prueba estadística.
3. el uso de la replicate función para que se ejecute el código de arriba un montón de veces (por lo general comienzan en alrededor de 100 veces para tener una idea de cuánto tiempo lleva y para obtener el derecho área general, hasta para 1,000 y, a veces, 10.000 o 100.000 para los valores finales que voy a utilizar). La proporción de veces que se rechazó la hipótesis nula es el poder.
4. rehacer el anterior para otro conjunto de condiciones.

Aquí es un simple ejemplo de 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 )

Answer 2

Además de la Nieve excelente ejemplo, yo creo que usted también puede hacer un poder de simulación por remuestreo de un conjunto de datos existente que tiene su efecto. No es un bootstrap, ya que no estás de muestreo con reemplazo de la misma n, pero la misma idea.

Así que he aquí un ejemplo: me encontré un poco de auto-experimento que se convirtió en un punto de partida positivo-estimación, sino porque era pequeño, no era casi estadísticamente significativa en el número ordinal de regresión logística. Con ese punto de estimar, de lo grande que un n necesito? Para varios posibles n, yo muchas veces genera un conjunto de datos y ejecutó el ordinal de regresión logística y vio cómo de pequeño es el p-valor fue de:

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 de (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, n=600 el poder fue del 32%. No es muy alentador.

(Si mi enfoque de simulación está mal, por favor que alguien me diga. Me voy un par de médicos de documentos de debate sobre el poder de la simulación para la planificación de ensayos clínicos, pero no estoy del todo seguro acerca de mi preciso de la aplicación.)

Answer 3

Me gustaría añadir una cosa a la Nieve de la respuesta (y esto se aplica a cualquier análisis del poder a través de la simulación) - prestar atención a si usted está buscando un 1 o 2 cola de la prueba. Programas populares como G*Potencia por defecto a 1 cola de la prueba, y si usted está tratando de ver si su simulaciones coinciden con ellos (siempre una buena idea cuando usted está aprendiendo cómo hacer esto), usted querrá comprobar que la primera.

Para hacer Nieve de ejecución de un 1 cola de prueba, me gustaría añadir un parámetro llamado "cola" a la función de las entradas y de poner algo como esto en la misma función:

 #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

El 1 de cola versión básicamente, se comprueba que el coeficiente es positivo y, a continuación, corta el p-valor en la mitad.