To pytanie jest odpowiedzią na odpowiedź udzieloną przez @Greg Snow na pytanie, które zadałem, dotyczące analizy mocy z regresją logistyczną i SAS Proc GLMPOWER
.
Jeśli projektuję eksperyment i przeanalizuję wyniki w silnej regresji logistycznej, jak mogę użyć symulacji (i tutaj ) do przeprowadzenia analizy mocy?
Oto prosty przykład, w którym istnieją dwie zmienne, pierwsza przyjmuje trzy możliwe wartości {0,03, 0,06, 0,09}, a druga jest fikcyjnym wskaźnikiem {0,1}. Dla każdego szacujemy współczynnik odpowiedzi dla każdej kombinacji (liczba respondentów / liczba osób wprowadzonych do obrotu). Ponadto chcemy mieć 3 razy więcej pierwszej kombinacji czynników niż innych (które można uznać za równe), ponieważ ta pierwsza kombinacja jest naszą wypróbowaną i prawdziwą wersją. Jest to konfiguracja podana w kursie SAS wspomnianym w powiązanym pytaniu.
Model, który zostanie wykorzystany do analizy wyników, będzie regresją logistyczną z głównymi efektami i interakcją (odpowiedź wynosi 0 lub 1).
mod <- glm(response ~ Var1 + Var2 + I(Var1*Var2))
Jak mogę zasymulować zestaw danych do użycia z tym modelem do przeprowadzenia analizy mocy?
Kiedy uruchomię to przez SAS Proc GLMPOWER
(użycie, STDDEV =0.05486016
które odpowiada sqrt(p(1-p))
gdzie p jest średnią ważoną pokazanych wskaźników odpowiedzi):
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;
Uwaga: GLMPOWER
będą używane tylko zmienne klasowe (nominalne), więc 3, 6, 9 powyżej są traktowane jako znaki i mogą być niskie, średnie i wysokie lub dowolne trzy pozostałe ciągi. Gdy przeprowadzana jest prawdziwa analiza, Var1 zostanie użyty jako wartość liczbowa (i uwzględnimy wielomianowy termin Var1 * Var1), aby uwzględnić dowolną krzywiznę.
Dane wyjściowe z SAS to
Widzimy więc, że potrzebujemy 762,112 jako naszej próbki (główny efekt Var2 jest najtrudniejszy do oszacowania) o mocy równej 0,80 i alfa równej 0,05. Przydzielilibyśmy je tak, aby 3 razy więcej było kombinacją wyjściową (tj. 0,375 * 762112), a reszta po prostu należała do pozostałych 5 kombinacji.
Odpowiedzi:
Czynności wstępne:
Jak omówiono w instrukcji G * Power , istnieje kilka różnych rodzajów analiz mocy, w zależności od tego, co chcesz rozwiązać. (To znaczy, , wielkość efektu , i moc istnieją w stosunku do siebie; określenie dowolnych trzech z nich pozwoli ci rozwiązać czwarty). E S αN ES α
Oprócz doskonałego posta @ GregSnow , kolejny naprawdę świetny przewodnik po analizach mocy opartych na symulacji na CV można znaleźć tutaj: Obliczanie mocy statystycznej . Podsumowując podstawowe pomysły:
To, czy znajdziesz znaczenie dla danej iteracji, może być rozumiane jako wynik próby Bernoulliego z prawdopodobieństwem (gdzie jest potęgą). Proporcja znaleziona w stosunku do iteracji pozwala nam zbliżyć się do prawdziwego . Aby uzyskać lepsze przybliżenie, możemy zwiększyć , chociaż spowoduje to również dłuższą symulację. p B p Bp p B p B
W języku R głównym sposobem generowania danych binarnych z danym prawdopodobieństwem „sukcesu” jest rbinom
rbinom(n=10, size=1, prob=p)
(prawdopodobnie będziesz chciał przypisać wynik do zmiennej do przechowywania)ifelse(runif(1)<=p, 1, 0)
pnorm()
i użyć tych wrbinom()
kodzie.Oświadczasz, że „uwzględnisz wielomianowy termin Var1 * Var1) w celu uwzględnienia dowolnej krzywizny”. Jest tu zamieszanie; warunki wielomianowe mogą pomóc nam uwzględnić krzywiznę, ale jest to termin interakcji - nie pomoże nam w ten sposób. Niemniej jednak wskaźniki odpowiedzi wymagają od nas uwzględnienia warunków kwadratowych i warunków interakcji w naszym modelu. W szczególności twój model będzie musiał zawierać: , i , poza podstawowymi terminami. v a r 1 ∗ v a r 2 v a r 1 2 ∗ v a r 2var12 var1∗var2 var12∗var2
Tak jak istnieją różne rodzaje wskaźników błędów typu I , gdy istnieje wiele hipotez (np wskaźnik błędu per-kontrast , poziom błędu familywise , i stopa błędów per-rodzina ), więc istnieją różne rodzaje zasilania * (na przykład do A jeden z góry określony efekt , dla dowolnego efektu i dla wszystkich efektów ). Możesz także szukać mocy do wykrycia określonej kombinacji efektów lub mocy jednoczesnego testu modelu jako całości. Domyślam się z twojego opisu twojego kodu SAS, że szuka tego drugiego. Jednak z twojego opisu twojej sytuacji zakładam, że chcesz przynajmniej wykryć efekty interakcji.
Aby zobaczyć inny sposób myślenia o kwestiach związanych z mocą, zobacz moją odpowiedź tutaj: Jak zgłosić ogólną precyzję w szacowaniu korelacji w kontekście uzasadnienia wielkości próby.
Prosta potęga post-hoc dla regresji logistycznej w R:
Powiedzmy, że podane przez ciebie wskaźniki odpowiedzi reprezentują prawdziwą sytuację na świecie i że wysłałeś 10 000 listów. Jaka jest moc wykrywania tych efektów? (Zauważ, że jestem znany z pisania „komicznie nieefektywnego” kodu, poniższe mają być łatwe do naśladowania, a nie zoptymalizowane pod kątem wydajności; w rzeczywistości są dość powolne).
Widzimy więc, że 10 000 liter tak naprawdę nie osiąga 80% mocy (jakiegokolwiek rodzaju) w celu wykrycia tych wskaźników odpowiedzi. (Nie jestem wystarczająco pewien, co robi kod SAS, aby móc wyjaśnić surową rozbieżność między tymi podejściami, ale kod ten jest koncepcyjnie prosty - jeśli jest wolny - i spędziłem trochę czasu na sprawdzaniu go, i myślę, że wyniki są rozsądne).
Moc a-priori oparta na symulacji dla regresji logistycznej:
Stąd pomysł polega na szukaniu możliwych , dopóki nie znajdziemy wartości, która daje pożądany poziom mocy, którą jesteś zainteresowany. Każda strategia wyszukiwania, którą możesz zakodować, aby z tym pracować, byłaby w porządku (w teoria). Biorąc pod uwagę , które będą wymagane do uchwycenia tak małych efektów, warto pomyśleć o tym, jak zrobić to bardziej efektywnie. Moje typowe podejście to po prostu brutalna siła, tj. Ocena każdego który mógłbym rozsądnie rozważyć. (Zauważ jednak, że zazwyczaj rozważałbym tylko niewielki zakres i zwykle pracuję z bardzo małymi wartościami - przynajmniej w porównaniu z tym.) N N NN N N N
Zamiast tego, moją strategią było ujęcie możliwych , aby zorientować się, jaki byłby zakres mocy. Tak więc wybrałem na 500 000 i ponownie uruchomiłem kod (inicjując to samo ziarno, ale uruchomienie trwało półtorej godziny). Oto wyniki: NN N
Widzimy z tego, że wielkość twoich efektów różni się znacznie, a zatem twoja zdolność do ich wykrywania jest różna. Na przykład efekt jest szczególnie trudny do wykrycia, stanowiąc jedynie 6% czasu, nawet przy pół miliona liter. Z drugiej strony, model jako całość był zawsze znacznie lepszy niż model zerowy. Inne możliwości są umieszczone pomiędzy. Mimo że większość „danych” jest wyrzucana podczas każdej iteracji, nadal możliwe jest spore badanie. Na przykład, moglibyśmy użyć macierzy do oceny korelacji między prawdopodobieństwami różnych zmiennych, które są znaczące.var12
significant
Na zakończenie powinienem zauważyć, że ze względu na złożoność i duże związane z twoją sytuacją, nie było to tak proste, jak podejrzewałem / twierdziłem w moim pierwszym komentarzu. Jednak na pewno możesz dowiedzieć się, jak to zrobić w ogóle i problemy związane z analizą mocy, z tego, co tu przedstawiłem. HTH.N
źródło
pwr
pakiet). Takie podejście daje ci możliwość znacznie wyraźniejszego myślenia (i / lub dopracowania swojego myślenia) o tym, czego oczekujesz, jak sobie z tym poradzisz itp. Uwaga: jednak potrzebujesz kwadratowych warunków lub czegoś takiego analogicznie, jeśli twoje ustawione prędkości są prawidłowe, b / c nie są liniowe, a sama interakcja nie pozwala uchwycić krzywoliniowych relacji.poly
zamiast pokazywania nowym użytkownikom R bardziej podatnej na błędy strategii kwadraturowania surowych wartości. Myślę, że pełny model powinien był zostaćglm(responses~ poly(var1, 2) * var2, family=binomial(link="logit")
. Byłby on zarówno mniej podatny na błędy statystyczne w interpretacji, jak i znacznie bardziej zwarty. Może nie być ważne w tym konkretnym przypadku, gdy patrzysz tylko na ogólne dopasowanie, ale może łatwo wprowadzić w błąd mniej wyrafinowanych użytkowników, którzy mogą ulec pokusie, aby spojrzeć na poszczególne warunki.=
itp. Ludzie będą bardziej zaznajomieni z kwadratowaniem zmiennych z podstawowej regresji klasy i mniej świadomi tego, copoly()
jest, jeśli nie są użytkownikami R.Odpowiedź @ Gunga jest świetna do zrozumienia. Oto podejście, które chciałbym zastosować:
Ta funkcja testuje ogólny efekt v2, modele można zmienić, aby spojrzeć na inne typy testów. Lubię pisać to jako funkcję, więc gdy chcę przetestować coś innego, mogę po prostu zmienić argumenty funkcji. Możesz także dodać pasek postępu lub użyć pakietu równoległego, aby przyspieszyć.
Tutaj zrobiłem 100 powtórzeń, zwykle zaczynam od tego poziomu, aby znaleźć przybliżoną wielkość próbki, a następnie zwiększam iteracje, gdy jestem w odpowiednim parku piłkarskim (nie musisz tracić czasu na 10 000 iteracji, gdy masz 20% mocy).
źródło
n
wierszami. Bardziej efektywne może być również wykonywanie wag w funkcji.