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

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.

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$$\alpha$

  • w swoim opisie chcesz znać odpowiedni aby uchwycić współczynniki odpowiedzi określone za pomocą i power = 80%. To jest moc a priori . α = $N$$\alpha=.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 tego$B$$N$
  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 $p$$p$$B$$p$$B$

• 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 2 ∗ v a r $var1^2$$var1*var2$$var1^2*var2$

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

  • * odniesienie: Maxwell, SE (2004). Trwałość słabych badań w badaniach psychologicznych: przyczyny, konsekwencje i środki zaradcze. Metody psychologiczne , 9 , 2 , s. 147-163.
  • twoje efekty są dość małe (nie mylić z niskimi wskaźnikami odpowiedzi), więc trudno będzie nam osiągnąć dobrą moc.
  • Zauważ, że chociaż wszystkie one brzmią dość podobnie, nie są bardzo takie same (np. Bardzo możliwe jest uzyskanie znaczącego modelu bez znaczących efektów - omówiono tutaj: Jak regresja może być znacząca, ale wszystkie predyktory mogą być inne niż znaczące? lub znaczące efekty, ale tam, gdzie model nie jest znaczący - omówiono tutaj: Istotność współczynników w regresji liniowej: znaczący test t vs nieistotna statystyka F ), co zostanie zilustrowane poniżej.

• 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 $N$$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: $N$$N$

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. $var1^2$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$

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