Analyse de puissance pour la régression logistique ordinale

12

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.

Peter Flom - Réintégrer Monica
la source

Réponses:

27

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 )
Greg Snow
la source
6
+1, il s'agit d'une approche universelle très robuste. Je l'ai souvent utilisé. Je voudrais suggérer une autre fonctionnalité: vous pouvez générer des données pour le maximum que vous considéreriez, puis ajuster le modèle pour les proportions de ces données en ajustant successivement le 1er n d'entre elles à intervalles réguliers jusqu'à (par exemple, n = 100 , 120, 140, 160, 180 et 200). Au lieu d'enregistrer une valeur p de chaque ensemble de données généré, vous pouvez enregistrer une ligne de valeurs p. La moyenne sur chaque colonne vous donne une idée rapide et sale de la façon dont la puissance change w / , et vous aide à affiner rapidement une valeur appropriée. NNN
gung - Rétablir Monica
2
@gung: votre commentaire a du sens, cela vous dérangerait-il d'ajouter vos codes afin que les personnes moins expérimentées en R puissent également en bénéficier? merci
1
Je regarde à nouveau ceci et j'ai quelques questions: 1) Pourquoi x uniforme est-il à 1:10? 2) Comment le généraliseriez-vous à plus d'une variable indépendante?
Peter Flom - Réintègre Monica
1
@PeterFlom, x devait être quelque chose, j'ai donc choisi (arbitrairement) de l'avoir uniforme entre 0 et 10, cela aurait aussi pu être normal, gamma, etc. Le mieux serait de choisir quelque chose qui est similaire à ce que nous attendons du réel x variables à ressembler. Pour utiliser plus d'une variable prédictive, générez-les indépendamment (ou à partir d'une normale multivariée, d'une copule, etc.), puis incluez-les toutes dans la pièce eta1, par exemple eta1 <- beta0 + beta1*x1 + beta2*x2 + beta3*x3.
Greg Snow
1
@ABC, ne pas répliquer ne vous donnerait qu'une seule décision, vous avez besoin de réplications pour déterminer la fréquence de rejet du test (la définition de la puissance). replicaten'est pas dans la fonction et ne modifie pas. La fonction renvoie la valeur de p (ce qui est en forme $ stats [5]) pour une itération, la réplication exécute la fonction 1000 fois (ou le nombre que vous spécifiez) et renvoie les 1000 valeurs de p, la meanfonction calcule ensuite la proportion de tests qui rejetteraient la valeur nulle à . α=0,05
Greg Snow
3

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é.

robin.datadrivers
la source
2

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

gwern
la source
0

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?

Karolina
la source
1
Je pense que vous obtiendriez une meilleure réponse si vous posiez une nouvelle question avec un lien vers cette question.
mdewey