Simulation logistischer Regressionskraftanalysen - entworfene Experimente

39

Diese Frage beantwortet @Greg Snow in Bezug auf eine Frage, die ich bezüglich der Leistungsanalyse mit logistischer Regression und SAS gestellt habe Proc GLMPOWER.

Wenn ich ein Experiment entwerfe und die Ergebnisse in einer faktoriellen logistischen Regression auswerten möchte, wie kann ich mithilfe der Simulation (und hier ) eine Leistungsanalyse durchführen?

Hier ist ein einfaches Beispiel mit zwei Variablen: Die erste nimmt drei mögliche Werte an {0,03, 0,06, 0,09} und die zweite ist ein Dummy-Indikator {0,1}. Für jede schätzen wir die Rücklaufquote für jede Kombination (Anzahl der Antwortenden / Anzahl der vermarkteten Personen). Darüber hinaus möchten wir, dass die erste Kombination von Faktoren dreimal so oft vorkommt wie die anderen (was als gleich angesehen werden kann), da diese erste Kombination unsere bewährte Version ist. Dies ist ein Setup, wie es in dem SAS-Kurs angegeben ist, der in der verknüpften Frage erwähnt wird.

Bildbeschreibung hier eingeben

Das Modell, das zur Analyse der Ergebnisse verwendet wird, ist eine logistische Regression mit Haupteffekten und Wechselwirkungen (Antwort ist 0 oder 1).

mod <- glm(response ~ Var1 + Var2 + I(Var1*Var2))

Wie kann ich einen Datensatz simulieren, der mit diesem Modell zur Durchführung einer Leistungsanalyse verwendet werden soll?

Wenn ich dies durch SAS ausgeführt Proc GLMPOWER(unter Verwendung von STDDEV =0.05486016 dem entspricht sqrt(p(1-p))den p der gewichtete Durchschnitt der gezeigten Ansprechraten ist):

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;

Hinweis: GLMPOWERNur Klassenvariablen (nominal) werden verwendet, sodass 3, 6, 9 oben als Zeichen behandelt werden und niedrig, mittel und hoch oder beliebige andere drei Zeichenfolgen sein können. Wenn die reelle Analyse durchgeführt wird, wird Var1 numerisch verwendet (und wir werden einen Polynomterm Var1 * Var1 einfügen), um jegliche Krümmung zu berücksichtigen.

Die Ausgabe von SAS ist

Bildbeschreibung hier eingeben

Wir sehen also, dass wir 762.112 als Stichprobengröße (Var2-Haupteffekt ist am schwierigsten abzuschätzen) mit einer Potenz von 0,80 und einem Alpha von 0,05 benötigen. Wir würden diese so zuordnen, dass die Basiskombination dreimal so groß ist (dh 0,375 * 762112) und der Rest nur zu gleichen Teilen in die anderen 5 Kombinationen fällt.

B_Miner
quelle
Dies ist in R einfach zu bewerkstelligen. 1. Frage: Stimmt es, dass 75% aller Fälle {var1 = .03, var2 = 0} & 25% für alle anderen Combos & nicht 3 Einheiten für jede 1 sein sollen? Einheit in jeder der anderen Combos (dh 37,5%)? 2. Frage, können Sie die Effekte angeben, die Sie erkennen möchten? Dh, was wäre die Log-Quote von 1 gegen 0? Wie sollten sich die Log-Erfolgsaussichten ändern, wenn var1 um 0,01 steigt? Glauben Sie, dass es zu einer Interaktion kommen könnte (wenn ja, wie groß ist die Interaktion?) (NB, diese Fragen können schwer zu beantworten sein. Eine Strategie besteht darin, den Anteil der
Einsen
1. Die Gewichtung von 3 für den Basisfall ist, dass es dreimal so viele Fälle gibt, in denen {var1 = 0,03, var2 = 0} ist. Die Ergebnisse von SAS (die besagen, dass wir 762.112 Stichprobengesamtgröße benötigen, um 80% der Hauptwirkung var2 = 0 zurückzuweisen, also die Gesamtstichprobengröße, die wir benötigen) würden diesem Basisfall 37,5% zugewiesen.
B_Miner
2. Nun, wir haben nur die Rücklaufquoten (das ist das erwartete Verhältnis der Anzahl der Erfolge zur Anzahl der Versuche). Wenn wir also 1.000 Briefe mit Var1 = 0,03 und Var2 = 0 senden, könnte dies einem Zinsangebot für ein Direktmailing-Angebot mit Kreditkarte von 0,03 (3%) und ohne Aufkleber auf dem Umschlag entsprechen (wobei Var2 = 1 bedeutet, dass dies der Fall ist) einen Aufkleber) erwarten wir 1000 * 0,0025 Antworten.
B_Miner
2. Forts.: Wir erwarten eine Interaktion - daher die Rücklaufquoten. Beachten Sie, dass für Var2 = 0 je nach Wert von Var1 eine andere Antwortrate gilt. Ich bin mir nicht sicher, wie ich diese in Log Odds und dann in einen simulierten Datensatz übersetzen soll.
B_Miner
Eine letzte Sache jedoch. Ich stelle fest, dass die Antwortraten für var1 linear sind, wenn var2 = 0 (dh .25%, .30%, .35%). Wollten Sie, dass dies ein linearer oder krummliniger Effekt ist? Sie sollten wissen, dass Wahrscheinlichkeiten für kleine Teilmengen ihres Bereichs ziemlich linear aussehen können, aber tatsächlich nicht linear sein können. Die logistische Regression wird linear in log odds, nicht Wahrscheinlichkeit (ich Sachen wie das in meiner Antwort diskutieren hier ).
gung - Wiedereinsetzung von Monica

Antworten:

43

Vorbereitungen:

  • NESα

    • Nα=.05
    • N
  • Neben dem hervorragenden Beitrag von @ GregSnow finden Sie hier einen weiteren wirklich tollen Leitfaden für simulationsbasierte Leistungsanalysen im CV: Berechnung der statistischen Leistung . Um die Grundideen zusammenzufassen:

    1. Ermitteln Sie den Effekt, den Sie erkennen möchten
    2. generiere N Daten aus dieser möglichen Welt
    3. Führen Sie die Analyse aus, die Sie für diese falschen Daten durchführen möchten
    4. Speichern Sie, ob die Ergebnisse für das von Ihnen gewählte Alpha "signifikant" sind
    5. BN
    6. N
  • ppBpB

  • In R ist der primäre Weg, um Binärdaten mit einer gegebenen Erfolgswahrscheinlichkeit zu erzeugen, das Binärsignal

    • Um beispielsweise die Anzahl der Erfolge aus 10 Bernoulli-Versuchen mit der Wahrscheinlichkeit p zu ermitteln, lautet der Code rbinom(n=10, size=1, prob=p)(Sie möchten das Ergebnis wahrscheinlich einer Variablen zur Speicherung zuweisen).
    • Sie können solche Daten auch weniger elegant erzeugen , indem Sie ? runif verwenden , z.ifelse(runif(1)<=p, 1, 0)
    • Wenn Sie glauben, dass die Ergebnisse durch eine latente Gaußsche Variable vermittelt werden, können Sie die latente Variable als Funktion Ihrer Kovariaten mit ? rnorm erzeugen und sie dann in Wahrscheinlichkeiten konvertieren pnorm()und diese in Ihrem rbinom()Code verwenden.
  • var12var1var2var12var2

  • Obwohl im Zusammenhang mit einer anderen Frage geschrieben, meine Antwort hier: Der Unterschied zwischen logit- und probit-Modellen enthält viele grundlegende Informationen zu diesen Modelltypen.
  • Ebenso wie es bei mehreren Hypothesen verschiedene Arten von Typ-I-Fehlerraten gibt (z. B. Fehlerrate pro Kontrast , familienweise Fehlerrate und familienweise Fehlerrate ), gibt es auch verschiedene Arten von Potenzen * (z. B. für a einen einzelnen vordefinierten Effekt für einen beliebigen Effekt und für alle Effekte ). Sie können auch nach der Fähigkeit suchen, eine bestimmte Kombination von Effekten zu erkennen, oder nach der Fähigkeit, das Modell als Ganzes gleichzeitig zu testen. Aus Ihrer Beschreibung Ihres SAS-Codes geht hervor, dass er nach letzterem sucht. Aus Ihrer Beschreibung Ihrer Situation gehe ich jedoch davon aus, dass Sie die Interaktionseffekte zumindest erfassen möchten.

  • Eine andere Art, über Fragen im Zusammenhang mit der Stromversorgung nachzudenken, finden Sie in meiner Antwort hier: Wie Sie die allgemeine Genauigkeit bei der Schätzung von Korrelationen im Kontext der Rechtfertigung der Stichprobengröße angeben.

Einfache Post-Hoc-Leistung für die logistische Regression in R:

Nehmen wir an, Ihre angegebenen Rücklaufquoten repräsentieren die wahre Situation auf der Welt und Sie haben 10.000 Briefe verschickt. Was ist die Kraft, um diese Effekte zu erkennen? (Beachten Sie, dass ich für das Schreiben von "komisch ineffizientem" Code berühmt bin. Das Folgende soll einfach zu befolgen sein, anstatt für die Effizienz optimiert zu sein. Tatsächlich ist es ziemlich langsam.)

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

Wir sehen also, dass 10.000 Briefe nicht wirklich 80% Leistung (jeglicher Art) erreichen, um diese Antwortraten zu erkennen. (Ich bin nicht sicher genug, was der SAS-Code tut, um die starke Diskrepanz zwischen diesen Ansätzen erklären zu können, aber dieser Code ist konzeptionell einfach - wenn auch langsam - und ich habe einige Zeit damit verbracht, ihn zu überprüfen, und ich denke, diese Ergebnisse sind vernünftig.)

Simulationsbasierte A-priori-Potenz für die logistische Regression:

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 - Wiedereinsetzung von Monica
quelle
Gung - WOW, vielen Dank für diese detaillierte und nachdenkliche Antwort! Wenn Sie meinen eigenen Code schreiben und mit Ihrem Code spielen, scheinen die quadratischen Terme das Problem zu sein - da bei einer viel kleineren Stichprobengröße eine Leistung von mindestens 80% erreicht wird, ohne dass dies im Modell berücksichtigt wird.
B_Miner
1
Das ist großartig, @B_Miner, genau das möchten Sie tun. Darüber hinaus ist der simulationsbasierte Ansatz meiner Meinung nach der analytischen Software überlegen, bei der nur eine Zahl ausgespuckt wird (R hat auch dieses pwrPaket). Dieser Ansatz gibt Ihnen die Möglichkeit, viel klarer darüber nachzudenken (und / oder Ihre Überlegungen zu verfeinern), was Sie erwarten, wie Sie damit umgehen würden Analog dazu sind Ihre gesetzten Raten, b / c, nicht linear, und die Interaktion allein lässt Sie keine krummlinigen Beziehungen erfassen.
gung - Wiedereinsetzung von Monica
Ich denke, Sie sollten die Verwendung der polyfehleranfälligeren Strategie des Quadrierens von Rohwerten demonstrieren, anstatt neuen Benutzern von R zu zeigen. Ich denke, das vollständige Modell hätte als gestellt werden sollen glm(responses~ poly(var1, 2) * var2, family=binomial(link="logit"). Es wäre sowohl weniger anfällig für statistische Interpretationsfehler als auch viel kompakter. In genau diesem Fall ist dies möglicherweise nicht wichtig, wenn Sie nur eine Gesamtanpassung betrachten, aber möglicherweise weniger erfahrene Benutzer irreführen, die möglicherweise versucht sind, sich einzelne Begriffe anzuschauen.
DW
@DWin, wenn ich R verwende, um die Dinge hier im Lebenslauf zu veranschaulichen, mache ich das auf eine sehr nicht-R-artige Weise. Die Idee ist, für diejenigen, die mit R nicht vertraut sind, so transparent wie möglich zu sein. ZB verwende ich nicht die vektorisierten Möglichkeiten, verwende Schleifen =usw. Die Leute sind mit dem Quadrieren von Variablen aus einer einfachen Regression besser vertraut Klasse, & weniger bewusst, was poly()ist, wenn sie nicht R-Benutzer sind.
gung - Reinstate Monica
17

@ Gungs Antwort ist großartig für das Verständnis. Hier ist der Ansatz, den ich verwenden würde:

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

Diese Funktion testet den Gesamteffekt von v2. Die Modelle können geändert werden, um andere Arten von Tests zu betrachten. Ich schreibe es gerne als Funktion, damit ich, wenn ich etwas anderes testen möchte, einfach die Funktionsargumente ändern kann. Sie können auch einen Fortschrittsbalken hinzufügen oder das Parallelpaket verwenden, um die Arbeit zu beschleunigen.

Hier habe ich gerade 100 Wiederholungen durchgeführt. Normalerweise beginne ich bei dieser Stufe, um die ungefähre Stichprobengröße zu ermitteln, und gehe dann die Iterationen hoch, wenn ich mich im richtigen Ballpark befinde (bei 20% Leistung muss ich keine Zeit mit 10.000 Iterationen verschwenden).

Greg Snow
quelle
Danke Greg! Ich habe mich über denselben Ansatz gewundert (wenn ich richtig verstehe, was Sie getan haben). Zur Bestätigung: Erstellen Sie einen Datensatz (in Kraft, aber mit Gewichten anstelle von Brute Force), der sehr groß ist und auf "mydat" basiert. , Anpassung einer logistischen Regression und anschließende Verwendung dieser Koeffizienten für die Stichprobe aus dem vorgeschlagenen Modell in der Simulation? Es scheint, dass dies ein allgemeiner Weg ist, um die Koeffizienten zu finden - dann ist es genau wie Ihre Antwort über die ordinale Regressionskraft, mit der ich verbunden bin.
B_Miner
Das ursprüngliche Modell verwendet Gewichte, um die zu verwendenden Koeffizienten zu ermitteln. In der Simulation wird jedoch ein Datenrahmen mit nZeilen erstellt. Es kann effizienter sein, Gewichte auch in der Funktion auszuführen.
Greg Snow
Ich habe Recht, dass Sie die Daten anfangs (um sehr gute Schätzungen zu erhalten) verwenden, um die verwendeten Koeffizienten zu ermitteln.
B_Miner
Ja, gut, das große ist so, dass das Verhältnis mal das Gewicht eine ganze Zahl ergibt.
Greg Snow
2
@B_Miner, ich plane einen Artikel, ich weiß nicht, ob es genug für ein vollständiges Buch gibt oder nicht. Aber danke für die Ermutigung.
Greg Snow