Analiza mocy dla regresji logistycznej porządkowej


Odpowiedzi:


27

Wolę robić analizy mocy poza podstawami poprzez symulację. W przypadku wstępnie zaplanowanych pakietów nigdy nie jestem całkiem pewien, jakie są założenia.

Symulacja mocy jest dość prosta (i niedroga) przy użyciu R.

  1. zdecyduj, jak Twoim zdaniem powinny wyglądać Twoje dane i jak je przeanalizujesz
  2. napisz funkcję lub zestaw wyrażeń, które będą symulować dane dla danej relacji i wielkości próbki i przeprowadzą analizę (funkcja jest lepsza, ponieważ możesz przekształcić wielkość próbki i parametry w argumenty, aby ułatwić wypróbowanie różnych wartości). Funkcja lub kod powinny zwracać wartość p lub inną statystykę testową.
  3. użyj tej replicatefunkcji, aby uruchomić kod z góry kilka razy (zwykle zaczynam od około 100 razy, aby sprawdzić, ile czasu to zajmuje i uzyskać odpowiedni obszar ogólny, a następnie zwiększyć go do 1000, a czasem 10 000 lub 100 000 dla ostateczne wartości, których użyję). Odsetek odrzuconych hipotez zerowych stanowi moc.
  4. powtórz powyższe dla innego zestawu warunków.

Oto prosty przykład z regresją porządkową:

library(rms)

tmpfun <- function(n, beta0, beta1, beta2) {
    x <- runif(n, 0, 10)
    eta1 <- beta0 + beta1*x
    eta2 <- eta1 + beta2
    p1 <- exp(eta1)/(1+exp(eta1))
    p2 <- exp(eta2)/(1+exp(eta2))
    tmp <- runif(n)
    y <- (tmp < p1) + (tmp < p2)
    fit <- lrm(y~x)
    fit$stats[5]
}

out <- replicate(1000, tmpfun(100, -1/2, 1/4, 1/4))
mean( out < 0.05 )

6
N.N.N.

2
@gung: twój komentarz ma sens, czy mógłbyś dodać swoje kody, aby mniej doświadczeni ludzie w R mogliby również z niego skorzystać? dzięki

1
Patrzę na to jeszcze raz i mam kilka pytań: 1) Dlaczego x uniform jest w skali 1:10? 2) Jak uogólniałbyś ją na więcej niż 1 niezależną zmienną?
Peter Flom - Przywróć Monikę

1
@PeterFlom, x musiał być czymś, więc wybrałem (arbitralnie), aby był jednolity między 0 a 10, mógł być również normalny, gamma itp. Najlepiej byłoby wybrać coś podobnego do tego, czego oczekujemy x zmienne do wyglądu. Aby użyć więcej niż 1 zmiennej predykcyjnej, wygeneruj je niezależnie (lub z wielowymiarowej wartości normalnej, kopuły itp.), A następnie po prostu dołącz je wszystkie do elementu eta1, np eta1 <- beta0 + beta1*x1 + beta2*x2 + beta3*x3.
Greg Snow

1
replicatemeanα=0,05

3

Dodam jeszcze jedną rzecz do odpowiedzi Snow'a (i dotyczy to każdej analizy mocy za pomocą symulacji) - zwróć uwagę na to, czy szukasz testu z 1 czy 2 ogonami. Popularne programy, takie jak G * Power domyślnie test 1-tailed, a jeśli próbujesz sprawdzić, czy Twoje symulacje pasują do nich (zawsze dobry pomysł, gdy uczysz się, jak to zrobić), powinieneś to sprawdzić najpierw.

Aby Snow uruchomił test jednostronny, dodałbym parametr o nazwie „ogon” do wejść funkcji i umieściłem coś takiego w samej funkcji:

 #two-tail test
  if (tail==2) fit$stats[5]

  #one-tail test
  if (tail==1){
    if (fit$coefficients[5]>0) {
          fit$stats[5]/2
    } else 1

Wersja 1-tailed zasadniczo sprawdza, czy współczynnik jest dodatni, a następnie zmniejsza wartość p o połowę.


2

Oprócz doskonałego przykładu Snow'a, wierzę, że możesz również przeprowadzić symulację mocy, ponownie próbkując z istniejącego zestawu danych, który ma na ciebie wpływ. Niezupełnie bootstrap, ponieważ nie próbujesz z tym samym n , ale ten sam pomysł.

Oto przykład: przeprowadziłem mały eksperyment własny, który przyniósł pozytywną ocenę punktową, ale ponieważ był mały, nie był prawie statystycznie istotny w regresji logistycznej porządkowej. Z tego punktu-kosztorysowej, jak duże n będę potrzebował? Dla różnych możliwych n , wiele razy wygenerowałem zestaw danych i przeprowadziłem regresję logistyczną porządkową i zobaczyłem, jak mała jest wartość p :

library(boot)
library(rms)
npt <- read.csv("http://www.gwern.net/docs/nootropics/2013-gwern-noopept.csv")
newNoopeptPower <- function(dt, indices) {
    d <- dt[sample(nrow(dt), n, replace=TRUE), ] # new dataset, possibly larger than the original
    lmodel <- lrm(MP ~ Noopept + Magtein, data = d)
    return(anova(lmodel)[7])
}
alpha <- 0.05
for (n in seq(from = 300, to = 600, by = 30)) {
   bs <- boot(data=npt, statistic=newNoopeptPower, R=10000, parallel="multicore", ncpus=4)
   print(c(n, sum(bs$t<=alpha)/length(bs$t)))
}

Z wyjściem (dla mnie):

[1] 300.0000   0.1823
[1] 330.0000   0.1925
[1] 360.0000   0.2083
[1] 390.0000   0.2143
[1] 420.0000   0.2318
[1] 450.0000   0.2462
[1] 480.000   0.258
[1] 510.0000   0.2825
[1] 540.0000   0.2855
[1] 570.0000   0.3184
[1] 600.0000   0.3175

W tym przypadku przy n = 600 moc wynosiła 32%. Niezbyt zachęcające.

(Jeśli moje podejście do symulacji jest niewłaściwe, proszę, powiedz mi. Idę z kilkoma artykułami medycznymi omawiającymi symulację mocy do planowania badań klinicznych, ale wcale nie jestem pewien co do mojej dokładnej realizacji).


0

Nawiązując do pierwszej symulacji (zaproponowanej przez Snow; /stats//a/22410/231675 ):

Nadal nie jestem pewien, jak powinna wyglądać symulacja z bardziej (konkretnie trzema) zmiennymi niezależnymi. Rozumiem, że powinienem „włączyć je wszystkie do kawałka eta1, np. Eta1 <- beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3 '' (jak wspomniano powyżej). Ale nie wiem, jak dostosować pozostałe parametry w funkcji. Czy ktoś mógłby mi w tym pomóc?


1
Myślę, że uzyskałbyś lepszą odpowiedź, gdybyś zadał nowe pytanie z linkiem do tego pytania.
mdewey
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.