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.
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.
Odpowiedzi:
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.
replicate
funkcji, 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.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 )
eta1 <- beta0 + beta1*x1 + beta2*x2 + beta3*x3
.
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?