Analyse de puissance pour la régression logistique ordinale

Je recherche un programme (en R ou SAS ou autonome, s'il est gratuit ou peu coûteux) qui fera une analyse de puissance pour la régression logistique ordinale.

Je préfère faire des analyses de puissance au-delà des bases par simulation. Avec les packages pré-planifiés, je ne suis jamais tout à fait sûr des hypothèses avancées.

La simulation de la puissance est assez simple (et abordable) en utilisant R.

1. décider à quoi devraient ressembler vos données et comment les analyser
2. écrire une fonction ou un ensemble d'expressions qui simuleront les données pour une relation et une taille d'échantillon données et effectuer l'analyse (une fonction est préférable dans la mesure où vous pouvez transformer la taille de l'échantillon et les paramètres en arguments pour faciliter l'essai de différentes valeurs). La fonction ou le code doit renvoyer la valeur p ou une autre statistique de test.
3. utiliser la replicatefonction pour exécuter le code à partir de plusieurs fois (je commence généralement à environ 100 fois pour avoir une idée du temps qu'il faut et pour obtenir la bonne zone générale, puis jusqu'à 1000 et parfois 10000 ou 100000 pour le valeurs finales que j'utiliserai). La proportion de fois où vous avez rejeté l'hypothèse nulle est la puissance.
4. refaire ce qui précède pour un autre ensemble de conditions.

Voici un exemple simple avec régression ordinale:

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 )

J'ajouterais une autre chose à la réponse de Snow (et cela s'applique à toute analyse de puissance via la simulation) - faites attention à si vous recherchez un test à 1 ou 2 queues. Des programmes populaires comme G * Power par défaut au test unilatéral, et si vous essayez de voir si vos simulations correspondent (toujours une bonne idée lorsque vous apprenez à le faire), vous voudrez d'abord vérifier cela.

Pour que Snow exécute un test unilatéral, j'ajouterais un paramètre appelé "queue" aux entrées de la fonction, et mettrais quelque chose comme ceci dans la fonction elle-même:

 #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 version unilatérale vérifie essentiellement que le coefficient est positif, puis réduit la valeur de p de moitié.

Outre l'excellent exemple de Snow, je pense que vous pouvez également faire une simulation de puissance en rééchantillonnant à partir d'un ensemble de données existant qui a votre effet. Pas tout à fait un bootstrap, car vous n'échantillonnez pas avec remplacement le même n , mais la même idée.

Voici donc un exemple: j'ai effectué une petite auto-expérience qui a donné une estimation ponctuelle positive, mais parce qu'elle était petite, n'était pas statistiquement significative dans la régression logistique ordinale. Avec ce point-estimation, la taille d' un n aurais - je besoin? Pour divers n possibles , j'ai à plusieurs reprises généré un ensemble de données et exécuté la régression logistique ordinale et vu à quel point la valeur p était petite :

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)))
}

Avec la sortie (pour moi):

[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

Dans ce cas, à n = 600, la puissance était de 32%. Pas très encourageant.

(Si mon approche de la simulation est erronée, s'il vous plaît, dites-le-moi. Je vais publier quelques articles médicaux sur la simulation de puissance pour planifier des essais cliniques, mais je ne suis pas du tout certain de ma mise en œuvre précise.)

En référence à la première simulation (suggérée par Snow; /stats//a/22410/231675 ):

Je ne sais toujours pas à quoi devrait ressembler la simulation avec plus (spécifiquement, trois) variables indépendantes. Je comprends que je devrais «les inclure tous dans la pièce eta1, par exemple eta1 <- beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3 '' (comme mentionné ci-dessus). Mais je ne sais pas comment régler le reste des paramètres de la fonction. Quelqu'un pourrait-il m'aider avec ça?