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.
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: GLMPOWER
Nur 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
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.
Antworten:
Vorbereitungen:
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:
In R ist der primäre Weg, um Binärdaten mit einer gegebenen Erfolgswahrscheinlichkeit zu erzeugen, das Binärsignal
rbinom(n=10, size=1, prob=p)
(Sie möchten das Ergebnis wahrscheinlich einer Variablen zur Speicherung zuweisen).ifelse(runif(1)<=p, 1, 0)
pnorm()
und diese in Ihremrbinom()
Code verwenden.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.)
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:
significant
quelle
pwr
Paket). 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.poly
fehleranfä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 sollenglm(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.=
usw. Die Leute sind mit dem Quadrieren von Variablen aus einer einfachen Regression besser vertraut Klasse, & weniger bewusst, waspoly()
ist, wenn sie nicht R-Benutzer sind.@ Gungs Antwort ist großartig für das Verständnis. Hier ist der Ansatz, den ich verwenden würde:
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).
quelle
n
Zeilen erstellt. Es kann effizienter sein, Gewichte auch in der Funktion auszuführen.