Cette question répond à une réponse donnée par @Greg Snow à une question que j’avais posée concernant l’analyse de puissance avec régression logistique et SAS Proc GLMPOWER
.
Si je conçois une expérience et que j'analyserai les résultats dans une régression logistique factorielle, comment puis-je utiliser la simulation (et ici ) pour effectuer une analyse de puissance?
Voici un exemple simple où il y a deux variables, la première prend trois valeurs possibles {0.03, 0.06, 0.09} et la seconde est un indicateur factice {0,1}. Pour chacun, nous estimons le taux de réponse pour chaque combinaison (nombre de répondants / nombre de personnes commercialisées). De plus, nous souhaitons avoir 3 fois plus de la première combinaison de facteurs que les autres (qui peuvent être considérés comme égaux) car cette première combinaison est notre version éprouvée. Ceci est une configuration comme celle donnée dans le cours SAS mentionné dans la question liée.
Le modèle qui sera utilisé pour analyser les résultats sera une régression logistique, avec des effets principaux et une interaction (la réponse est 0 ou 1).
mod <- glm(response ~ Var1 + Var2 + I(Var1*Var2))
Comment puis-je simuler un ensemble de données à utiliser avec ce modèle pour effectuer une analyse de puissance?
Lorsque je lance ceci via SAS Proc GLMPOWER
(en utilisant STDDEV =0.05486016
ce qui correspond à sqrt(p(1-p))
où p est la moyenne pondérée des taux de réponse indiqués):
data exemplar;
input Var1 $ Var2 $ response weight;
datalines;
3 0 0.0025 3
3 1 0.00395 1
6 0 0.003 1
6 1 0.0042 1
9 0 0.0035 1
9 1 0.002 1;
run;
proc glmpower data=exemplar;
weight weight;
class Var1 Var2;
model response = Var1 | Var2;
power
power=0.8
ntotal=.
stddev=0.05486016;
run;
Remarque: GLMPOWER
seules les variables de classe (nominales) seront utilisées. Ainsi, 3, 6 et 9 ci-dessus sont traités comme des caractères et peuvent être constitués des chaînes basse, moyenne et haute, ou de trois autres chaînes. Lorsque la véritable analyse est effectuée, Var1 sera utilisé un numérique (et nous inclurons un terme polynomial Var1 * Var1) pour rendre compte de toute courbure.
La sortie de SAS est
Nous voyons donc que nous avons besoin de 762 112 comme la taille de notre échantillon (l’effet principal Var2 est le plus difficile à estimer) avec une puissance égale à 0,80 et une valeur alpha égale à 0,05. Nous attribuerions ces ressources à trois fois plus que la combinaison de base (soit 0,375 * 762112) et le reste se répartit également entre les 5 autres combinaisons.
Réponses:
Préliminaires:
En plus de l' excellent article de @ GregSnow, vous trouverez un autre excellent guide pour l'analyse de la puissance sur CV basée sur la simulation: Calcul de la puissance statistique . Pour résumer les idées de base:
En R, le premier moyen de générer des données binaires avec une probabilité de «succès» donnée est ? Rbinom
rbinom(n=10, size=1, prob=p)
(vous voudrez probablement affecter le résultat à une variable pour le stockage).ifelse(runif(1)<=p, 1, 0)
pnorm()
celles de votrerbinom()
code et les utiliser .Tout comme il existe différents types de taux d'erreurs I quand il y a plusieurs hypothèses (par exemple, le taux d'erreur par contraste , le taux d'erreur de l' , et par famille taux d'erreur ), alors qu'il existe différents types de pouvoir * (par exemple, pour un effet unique prédéterminé , pour tout effet , et pour tous les effets ). Vous pouvez également rechercher le pouvoir de détecter une combinaison d'effets spécifique ou le pouvoir d'effectuer un test simultané du modèle dans son ensemble. Je suppose que d'après votre description de votre code SAS, il recherche ce dernier. Cependant, d'après votre description de votre situation, je suppose que vous souhaitez détecter les effets d'interaction au minimum.
Pour une manière différente de penser aux problèmes liés au pouvoir, voir ma réponse ici: Comment rendre compte de la précision générale dans l'estimation des corrélations dans un contexte de justification de la taille de l'échantillon.
Puissance post-hoc simple pour la régression logistique dans R:
Supposons que vos taux de réponse supposés représentent la situation réelle dans le monde et que vous avez envoyé 10 000 lettres. Quel est le pouvoir de détecter ces effets? (Notez que je suis célèbre pour l'écriture de code "comiquement inefficace", ce qui suit est destiné à être facile à suivre plutôt qu'optimisé pour plus d'efficacité; en fait, c'est assez lent.)
Nous voyons donc que 10 000 lettres n'atteignent pas 80% de la puissance (de quelque sorte que ce soit) pour détecter ces taux de réponse. (Je ne suis pas assez sûr de ce que fait le code SAS pour pouvoir expliquer la contradiction flagrante entre ces approches, mais ce code est conceptuellement simple - même s’il est lent - et j’ai passé un certain temps à le vérifier, et je pense que ces les résultats sont raisonnables.)
Puissance a priori basée sur la simulation pour la régression logistique:
significant
la source
pwr
paquet). Cette approche vous donne l’occasion de penser beaucoup plus clairement (et / ou d’affiner votre réflexion) sur ce que vous attendez, sur la façon dont vous traitez cela, etc. Remarquez cependant que vous avez besoin des termes quadratiques, ou de quelque chose de différent. de façon analogue, si vos tarifs proposés sont corrects, car ils ne sont pas linéaires, et l’interaction seule ne vous permet pas de capturer des relations curvilignes.poly
plutôt que de montrer aux nouveaux utilisateurs de R la stratégie plus prédisposée aux erreurs de quadrature des valeurs brutes. Je pense que le modèle complet aurait dû être posé commeglm(responses~ poly(var1, 2) * var2, family=binomial(link="logit")
. Ce serait à la fois moins sujet aux erreurs statistiques d’interprétation et beaucoup plus compact. Cela pourrait ne pas être important dans ce cas précis lorsque vous ne tenez compte que de l’ajustement global, mais que vous risquez d’induire en erreur des utilisateurs moins avertis qui pourraient être tentés d’examiner des termes individuels.=
, etc. Les gens seront plus familiarisés avec les variables de quadrature à partir d'une régression de base. classe, et moins au courant de ce quipoly()
est, s’ils ne sont pas des utilisateurs de R.La réponse de @ Gung est excellente pour la compréhension. Voici l'approche que j'utiliserais:
Cette fonction teste l'effet global de v2, les modèles peuvent être modifiés pour examiner d'autres types de tests. J'aime l'écrire en tant que fonction pour pouvoir modifier les arguments de la fonction lorsque je souhaite tester quelque chose de différent. Vous pouvez également ajouter une barre de progression ou utiliser le package parallèle pour accélérer les choses.
Ici, je viens de faire 100 répétitions, je commence habituellement autour de ce niveau pour trouver la taille approximative de l’échantillon, puis monte les itérations lorsque je suis dans le bon stade (inutile de perdre du temps sur 10 000 itérations lorsque vous avez 20% de puissance).
la source
n
lignes. Il serait peut-être plus efficace de faire de la pondération dans la fonction également.