Analiza mocy dla regresji logistycznej porządkowej

Szukam programu (w wersji R lub SAS lub autonomicznej, jeśli jest darmowy lub tani), który wykona analizę mocy dla regresji logistycznej porządkowej.

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 )

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

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

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?