Simulation d'analyse de puissance de régression logistique - expériences conçues

39

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.

entrez la description de l'image ici

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: GLMPOWERseules 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

entrez la description de l'image ici

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.

B_Miner
la source
Ceci est facile à faire dans R. 1ère question: ai-je raison de dire que vous voulez que 75% de tous les cas soient {var1 = .03, var2 = 0} & 25% pour tous les autres combos, et non 3 unités pour chaque 1 unité dans chacun des autres combos (c.-à-d. 37,5%)? 2ème question, pouvez-vous préciser les effets que vous souhaitez détecter? C'est-à-dire, quelle serait la cote du journal de 1 contre 0? Comment les chances de réussite du journal devraient-elles changer si la var1 augmente de 0,01? Pensez-vous qu'il pourrait y avoir une interaction (si oui, quelle est sa taille)? (NB, il peut être difficile de répondre à ces questions, une stratégie consiste à spécifier la proportion de réponses que vous pensez être dans chaque combo.)
gung - Reinstate Monica
1er: Le poids de 3 pour le cas de base est qu'il y a 3 fois plus de cas où {var1 = 0.03, var2 = 0}. Ainsi, les résultats de SAS (qui dit que nous avons besoin de 762 112 tailles d'échantillon totales pour pouvoir rejeter l'effet principal à var2 = 0, soit la taille d'échantillon totale dont nous avons besoin) se verraient attribuer 37,5% à ce cas de base.
B_Miner
2ème: Tout ce que nous avons, c'est le taux de réponse (qui est le ratio attendu du nombre de succès sur le nombre d'essais). Donc, si nous envoyons 1 000 lettres avec Var1 = 0,03 et Var2 = 0, cela pourrait correspondre à une offre de taux d'intérêt sur une offre de publipostage par carte de crédit de 0,03 (3%) et l'absence d'autocollant sur l'enveloppe (où Var2 = 1 signifie qu'il y a un autocollant), nous attendons 1000 * 0,0025 réponses.
B_Miner
2ème cont: Nous attendons une interaction - d'où les taux de réponse. Notez qu'il existe un taux de réponse différent pour Var2 = 0 en fonction de la valeur de Var1. Je ne suis pas sûr de savoir comment traduire ces informations pour consigner les cotes, puis les données simulées.
B_Miner
Une dernière chose, cependant. Je remarque que les taux de réponse sont linéaires pour var1 lorsque var2 = 0 (soit 0,25%, 0,30%, 0,35%). Aviez-vous l'intention que cet effet soit linéaire ou curviligne? Vous devez savoir que les probabilités peuvent sembler assez linéaires pour de petits sous-ensembles de leur plage, mais ne peuvent en réalité être linéaires. La régression logistique est linéaire en cotes de log, pas en probabilité (je discute de ce genre de choses dans ma réponse ici ).
gung - Réintégrer Monica

Réponses:

43

Préliminaires:

  • NESα

    • Nα=.05
    • N
  • 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:

    1. comprendre l'effet que vous voulez pouvoir détecter
    2. générer N données de ce monde possible
    3. lancez l'analyse que vous avez l'intention de faire sur ces fausses données
    4. stocker si les résultats sont «significatifs» selon l'alpha choisi
    5. BN
    6. N
  • ppBpB

  • En R, le premier moyen de générer des données binaires avec une probabilité de «succès» donnée est ? Rbinom

    • Par exemple, pour obtenir le nombre de succès sur 10 essais de Bernoulli avec probabilité p, le code serait rbinom(n=10, size=1, prob=p)(vous voudrez probablement affecter le résultat à une variable pour le stockage).
    • vous pouvez également générer de telles données de manière moins élégante en utilisant ? runif , par exemple,ifelse(runif(1)<=p, 1, 0)
    • Si vous pensez que les résultats sont médiés par une variable gaussienne latente, vous pouvez générer la variable latente en fonction de vos covariables avec normrnorm , puis les convertir en probabilités avec pnorm()celles de votre rbinom()code et les utiliser .
  • var12var1var2var12var2

  • Bien que rédigé dans le contexte d’une autre question, ma réponse est la suivante: La différence entre les modèles logit et probit contient de nombreuses informations de base sur ces types de modèles.
  • 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.)

set.seed(1)

repetitions = 1000
N = 10000
n = N/8
var1  = c(   .03,    .03,    .03,    .03,    .06,    .06,    .09,   .09)
var2  = c(     0,      0,      0,      1,      0,      1,      0,     1)
rates = c(0.0025, 0.0025, 0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002)

var1    = rep(var1, times=n)
var2    = rep(var2, times=n)
var12   = var1**2
var1x2  = var1 *var2
var12x2 = var12*var2

significant = matrix(nrow=repetitions, ncol=7)

startT = proc.time()[3]
for(i in 1:repetitions){
  responses          = rbinom(n=N, size=1, prob=rates)
  model              = glm(responses~var1+var2+var12+var1x2+var12x2, 
                           family=binomial(link="logit"))
  significant[i,1:5] = (summary(model)$coefficients[2:6,4]<.05)
  significant[i,6]   = sum(significant[i,1:5])
  modelDev           = model$null.deviance-model$deviance
  significant[i,7]   = (1-pchisq(modelDev, 5))<.05
}
endT = proc.time()[3]
endT-startT

sum(significant[,1])/repetitions      # pre-specified effect power for var1
[1] 0.042
sum(significant[,2])/repetitions      # pre-specified effect power for var2
[1] 0.017
sum(significant[,3])/repetitions      # pre-specified effect power for var12
[1] 0.035
sum(significant[,4])/repetitions      # pre-specified effect power for var1X2
[1] 0.019
sum(significant[,5])/repetitions      # pre-specified effect power for var12X2
[1] 0.022
sum(significant[,7])/repetitions      # power for likelihood ratio test of model
[1] 0.168
sum(significant[,6]==5)/repetitions   # all effects power
[1] 0.001
sum(significant[,6]>0)/repetitions    # any effect power
[1] 0.065
sum(significant[,4]&significant[,5])/repetitions   # power for interaction terms
[1] 0.017

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:

NNNN

NN

sum(significant[,1])/repetitions      # pre-specified effect power for var1
[1] 0.115
sum(significant[,2])/repetitions      # pre-specified effect power for var2
[1] 0.091
sum(significant[,3])/repetitions      # pre-specified effect power for var12
[1] 0.059
sum(significant[,4])/repetitions      # pre-specified effect power for var1X2
[1] 0.606
sum(significant[,5])/repetitions      # pre-specified effect power for var12X2
[1] 0.913
sum(significant[,7])/repetitions      # power for likelihood ratio test of model
[1] 1
sum(significant[,6]==5)/repetitions   # all effects power
[1] 0.005
sum(significant[,6]>0)/repetitions    # any effect power
[1] 0.96
sum(significant[,4]&significant[,5])/repetitions   # power for interaction terms
[1] 0.606

var12significant

N

gung - Rétablir Monica
la source
Gung - WOW, merci beaucoup pour cette réponse détaillée et réfléchie! En écrivant le mien et en jouant avec votre code, les termes du second degré semblent être un problème, car au moins 80% de la puissance est atteinte avec une taille d'échantillon beaucoup plus petite, sans en tenir compte dans le modèle.
B_Miner
1
C'est génial, @B_Miner, c'est le genre de chose que vous voulez faire. De plus, c’est la raison pour laquelle je pense que l’approche basée sur la simulation est supérieure au logiciel analytique qui crache un nombre (R a cela aussi, le pwrpaquet). 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.
gung - Réintégrer Monica
Je pense que vous devriez démontrer l’utilisation de polyplutô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é comme glm(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.
DWin
@DWin, lorsque j'utilise R pour illustrer des choses ici sur CV, je le fais d'une manière très différente de R. L'idée est d'être aussi transparent que possible pour ceux qui ne sont pas familiers avec R. Par exemple, je n'utilise pas les possibilités vectorisées, j'utilise des boucles =, 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 qui poly()est, s’ils ne sont pas des utilisateurs de R.
gung - Rétablir Monica
17

La réponse de @ Gung est excellente pour la compréhension. Voici l'approche que j'utiliserais:

mydat <- data.frame( v1 = rep( c(3,6,9), each=2 ),
    v2 = rep( 0:1, 3 ), 
    resp=c(0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002) )

fit0 <- glm( resp ~ poly(v1, 2, raw=TRUE)*v2, data=mydat,
    weight=rep(100000,6), family=binomial)
b0 <- coef(fit0)


simfunc <- function( beta=b0, n=10000 ) {
    w <- sample(1:6, n, replace=TRUE, prob=c(3, rep(1,5)))
    mydat2 <- mydat[w, 1:2]
    eta <- with(mydat2,  cbind( 1, v1, 
                v1^2, v2,
                v1*v2,
                v1^2*v2 ) %*% beta )
    p <- exp(eta)/(1+exp(eta))
    mydat2$resp <- rbinom(n, 1, p)

    fit1 <- glm( resp ~ poly(v1, 2)*v2, data=mydat2,
        family=binomial)
    fit2 <- update(fit1, .~ poly(v1,2) )
    anova(fit1,fit2, test='Chisq')[2,5]
}

out <- replicate(100, simfunc(b0, 10000))
mean( out <= 0.05 )
hist(out)
abline(v=0.05, col='lightgrey')

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

Greg Snow
la source
Merci Greg! Je m'interrogeais sur cette même approche (si je comprends bien ce que vous avez fait). Pour confirmer: créez-vous un ensemble de données (en fait, mais avec des pondérations au lieu de la force brute créant des enregistrements individuels des valeurs de Var1 et Var2, puis des 1 et des 0 pour la réponse) très volumineux basé sur "mydat" , ajustant une régression logistique et utilisant ensuite ces coefficients pour échantillonner à partir du modèle proposé dans la simulation? Il semble que ce soit un moyen général d’arriver aux coefficients - c’est alors tout comme votre réponse au sujet du pouvoir de régression ordinale auquel je suis lié.
B_Miner
Le modèle initial utilise des poids pour obtenir les coefficients à utiliser, mais dans la simulation, il crée un bloc de données avec des nlignes. Il serait peut-être plus efficace de faire de la pondération dans la fonction également.
Greg Snow
J'ai raison de dire que vous utilisez initialement les données (ce qui vous permet d'obtenir de très bonnes estimations) pour obtenir les coefficients utilisés?
B_Miner
Oui, eh bien le gros est tel que la proportion multipliée par le poids donne un entier.
Greg Snow
2
@B_Miner, je prévois un article, je ne sais pas s'il y en a assez pour un livre complet ou pas. Mais merci pour les encouragements.
Greg Snow