Symulacja analizy mocy regresji logistycznej - zaprojektowane eksperymenty

39

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.

wprowadź opis zdjęcia tutaj

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: GLMPOWERbę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

wprowadź opis zdjęcia tutaj

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.

B_Miner
źródło
Łatwo to zrobić w R. Pierwsze pytanie: czy mam rację, że chcesz 75% wszystkich przypadków to {var1 = .03, var2 = 0} i 25% dla wszystkich innych kombinacji, a nie 3 jednostki na każde 1 jednostka w każdym z pozostałych zestawów (37,5%)? Drugie pytanie, czy możesz określić efekty, które chcesz wykryć? To znaczy, jakie byłyby dzienne szanse 1 na 0? Jak szanse na sukces dziennika powinny się zmienić, jeśli var1 wzrośnie o 0,01? Czy uważasz, że może wystąpić interakcja (jeśli tak, to jak duża)? (Uwaga, te Q mogą być trudne do odpowiedzi, 1 strategia polega na określeniu proporcji 1, które Twoim zdaniem mogą znajdować się w każdej kombinacji).
gung - Przywróć Monikę
Po pierwsze: waga 3 dla przypadku wyjściowego jest taka, że ​​istnieje 3 razy więcej przypadków, w których {var1 = 0,03, var2 = 0}. Tak więc wyniki z SAS (który mówi, że potrzebujemy 762,112 całkowitej wielkości próbki, aby mieć 80% mocy odrzucenia efektu głównego var2 = 0, więc to jest potrzebna całkowita wielkość próby), przydzielono by 37,5% do tego podstawowego scenariusza.
B_Miner
Po drugie: cóż, wszystko, co mamy, to wskaźniki odpowiedzi (czyli oczekiwany stosunek liczby sukcesów do liczby prób). Tak więc, jeśli wyślemy 1000 listów z Var1 = 0,03 i Var2 = 0, co może odpowiadać oprocentowaniu oferty pocztowej karty kredytowej w wysokości 0,03 (3%) i braku naklejki na kopercie (gdzie Var2 = 1 oznacza, że ​​istnieje naklejka), oczekujemy 1000 * 0,0025 odpowiedzi.
B_Miner
2. cd: Oczekujemy interakcji - stąd wskaźniki odpowiedzi. Zauważ, że istnieje inny wskaźnik odpowiedzi dla Var2 = 0 w zależności od wartości Var1. Nie jestem pewien, jak je przetłumaczyć, aby zapisać szanse, a następnie symulowany zestaw danych.
B_Miner
Ale ostatnia rzecz. Zauważam, że współczynniki odpowiedzi są liniowe dla var1, gdy var2 = 0 (tj. 0,25%, 0,30%, 0,35%). Czy chciałeś, aby był to efekt liniowy czy krzywoliniowy? Powinieneś wiedzieć, że prawdopodobieństwa mogą wyglądać dość liniowo dla małych podzbiorów ich zakresu, ale w rzeczywistości nie mogą być liniowe. Regresja logistyczna jest liniowa w logarytmicznych szansach, a nie prawdopodobieństwie (omawiam takie rzeczy w mojej odpowiedzi tutaj ).
gung - Przywróć Monikę

Odpowiedzi:

43

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 αNESα

    • w swoim opisie chcesz znać odpowiedni aby uchwycić współczynniki odpowiedzi określone za pomocą i power = 80%. To jest moc a priori . α = 0,05Nα=.05
    • możemy zacząć od mocy post-hoc (określić moc na podstawie , odpowiedzi i alfa), ponieważ jest to koncepcyjnie prostsze, a następnie przejść w góręN
  • 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:

    1. dowiedzieć się, jaki efekt chcesz wykryć
    2. generuj N danych z tego możliwego świata
    3. uruchom analizę, którą zamierzasz przeprowadzić na tych fałszywych danych
    4. zapisz, czy wyniki są „znaczące” zgodnie z wybraną przez ciebie alfą
    5. powtórz wiele ( ) razy i wykorzystaj% „znaczący” jako oszacowanie mocy (post-hoc) dla tegoNBN
    6. aby określić moc a priori, wyszukaj możliwe wartości , aby znaleźć wartość, która daje pożądaną moc N
  • 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 BppBpB

  • W języku R głównym sposobem generowania danych binarnych z danym prawdopodobieństwem „sukcesu” jest rbinom

    • Na przykład, aby uzyskać liczbę sukcesów z 10 prób Bernoulliego z prawdopodobieństwem p, kod byłby rbinom(n=10, size=1, prob=p)(prawdopodobnie będziesz chciał przypisać wynik do zmiennej do przechowywania)
    • możesz również generować takie dane mniej elegancko, używając ? runif , np.ifelse(runif(1)<=p, 1, 0)
    • jeśli uważasz, że wyniki są zapośredniczone przez ukrytą zmienną Gaussa, możesz wygenerować zmienną ukrytą jako funkcję twoich zmiennych towarzyszących za pomocą ? rnorm , a następnie przekształcić je w prawdopodobieństwa pnorm()i użyć tych w rbinom()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 2v a r 2var12var1var2var12var2

  • Chociaż napisane w kontekście innego pytania, moja odpowiedź tutaj: Różnica między modelami logit i probit zawiera wiele podstawowych informacji o tych typach modeli.
  • 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).

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

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 NNNNN

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

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

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

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

gung - Przywróć Monikę
źródło
Gung - WOW, dziękuję bardzo za tak szczegółową i przemyślaną odpowiedź! Pisząc mój własny i grając z twoim kodem, wydaje się, że problemem są kwadratowe terminy - ponieważ co najmniej 80% mocy jest uzyskiwane przy znacznie mniejszym rozmiarze próbki bez uwzględnienia tego w modelu.
B_Miner
1
To świetnie, @B_Miner, to jest coś, co chcesz robić. Co więcej, jest to powód, dla którego uważam, że podejście oparte na symulacji jest lepsze od oprogramowania analitycznego, które po prostu wyrzuca liczbę (R ma to także, pwrpakiet). 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.
gung - Przywróć Monikę
Myślę, że powinieneś demonstrować użycie polyzamiast 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.
DWin
@DWin, kiedy używam R do zilustrowania rzeczy tutaj na CV, robię to w bardzo nie-R sposób. Chodzi o to, aby być jak najbardziej przejrzystym dla tych, którzy nie są zaznajomieni z R. Na przykład, nie używam wektoryzowanych możliwości, używam pętli =itp. Ludzie będą bardziej zaznajomieni z kwadratowaniem zmiennych z podstawowej regresji klasy i mniej świadomi tego, co poly()jest, jeśli nie są użytkownikami R.
Gung - Przywróć Monikę
17

Odpowiedź @ Gunga jest świetna do zrozumienia. Oto podejście, które chciałbym zastosować:

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

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

Greg Snow
źródło
Dzięki Greg! Zastanawiałem się nad tym samym podejściem (jeśli dobrze rozumiem, co zrobiłeś). Aby potwierdzić: czy tworzysz zestaw danych (w rzeczywistości, ale robisz to z wagami zamiast brutalnej siły, tworząc indywidualne rekordy wartości Var1 i Var2, a następnie 1 i 0 dla odpowiedzi), który jest bardzo duży na podstawie „mydat” , dopasowanie regresji logistycznej, a następnie wykorzystanie tych współczynników do pobrania próbki z proponowanego modelu w symulacji? Wygląda na to, że jest to ogólny sposób na wyliczenie współczynników - to jest tak jak twoja odpowiedź na temat regresji porządkowej, z którą się połączyłem.
B_Miner
Początkowy model wykorzystuje wagi, aby uzyskać współczynniki do użycia, ale w symulacji tworzy ramkę danych z nwierszami. Bardziej efektywne może być również wykonywanie wag w funkcji.
Greg Snow,
Mam rację, że początkowo używasz danych (co czyni je dużymi, aby uzyskać bardzo dobre oszacowania) w celu uzyskania używanych współczynników?
B_Miner
Tak, cóż, duża jest taka, że ​​proporcja razy waga daje liczbę całkowitą.
Greg Snow,
2
@B_Miner, planuję artykuł, nie wiem, czy wystarczy na pełną książkę, czy nie. Ale dziękuję za zachętę.
Greg Snow