Analisis daya untuk regresi logistik ordinal

Saya mencari program (dalam R atau SAS atau mandiri, jika gratis atau murah) yang akan melakukan analisis daya untuk regresi logistik ordinal.

Saya lebih suka melakukan analisis kekuatan di luar dasar-dasar dengan simulasi. Dengan paket pra-paket, saya tidak pernah yakin asumsi apa yang dibuat.

Simulasi untuk daya cukup mudah (dan terjangkau) menggunakan R.

1. putuskan seperti apa data Anda seharusnya terlihat dan bagaimana Anda akan menganalisisnya
2. tulis fungsi atau serangkaian ekspresi yang akan mensimulasikan data untuk hubungan dan ukuran sampel tertentu dan melakukan analisis (fungsi lebih disukai karena Anda dapat membuat ukuran dan parameter sampel menjadi argumen untuk membuatnya lebih mudah untuk mencoba nilai yang berbeda). Fungsi atau kode harus mengembalikan nilai p atau statistik pengujian lainnya.
3. gunakan replicatefungsi untuk menjalankan kode dari atas beberapa kali (saya biasanya mulai sekitar 100 kali untuk merasakan berapa lama dan untuk mendapatkan area umum yang tepat, kemudian naik menjadi 1.000 dan kadang-kadang 10.000 atau 100.000 untuk nilai akhir yang akan saya gunakan). Proporsi kali Anda menolak hipotesis nol adalah kekuatan.
4. ulangi di atas untuk serangkaian kondisi lain.

Berikut adalah contoh sederhana dengan regresi ordinal:

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 )

Saya akan menambahkan satu hal lain untuk jawaban Snow (dan ini berlaku untuk analisis daya apa pun melalui simulasi) - perhatikan apakah Anda mencari tes 1 atau 2 tailed. Program-program populer seperti G * Power default ke tes 1-tailed, dan jika Anda mencoba untuk melihat apakah simulasi Anda cocok dengan mereka (selalu merupakan ide bagus ketika Anda belajar bagaimana melakukan ini), Anda akan ingin memeriksanya terlebih dahulu.

Untuk membuat Snow menjalankan tes 1-tailed, saya akan menambahkan parameter yang disebut "tail" ke input fungsi, dan meletakkan sesuatu seperti ini di fungsi itu sendiri:

 #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

Versi 1-tailed pada dasarnya memeriksa untuk melihat bahwa koefisiennya positif, dan kemudian memotong nilai-p menjadi dua.

Selain contoh bagus Snow, saya yakin Anda juga dapat melakukan simulasi daya dengan melakukan resampling dari dataset yang sudah ada yang memiliki efek Anda. Tidak cukup bootstrap, karena Anda tidak sampel-dengan-pengganti yang sama n , tapi ide yang sama.

Jadi inilah contohnya: Saya menjalankan sedikit eksperimen sendiri yang menghasilkan estimasi titik positif tetapi karena jumlahnya sedikit, hampir tidak signifikan secara statistik dalam regresi logistik ordinal. Dengan estimasi titik itu, seberapa besar dan n yang saya butuhkan? Untuk berbagai kemungkinan n , saya berkali-kali membuat dataset & menjalankan regresi logistik ordinal & melihat seberapa kecil p- value-nya:

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

Dengan output (untuk saya):

[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

Dalam hal ini, pada n = 600 daya adalah 32%. Tidak terlalu membesarkan hati.

(Jika pendekatan simulasi saya salah, tolong seseorang memberi tahu saya. Saya akan pergi beberapa makalah medis membahas simulasi daya untuk merencanakan uji klinis, tapi saya sama sekali tidak yakin tentang implementasi saya yang tepat.)

Mengacu pada simulasi pertama (disarankan oleh Snow; /stats//a/22410/231675 ):

Saya masih tidak yakin bagaimana simulasi akan terlihat seperti dengan lebih banyak (khususnya, tiga) variabel independen. Saya mengerti bahwa saya harus 'memasukkan semuanya ke dalam eta1 piece, mis. Eta1 <- beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3' '(seperti yang disebutkan di atas). Tapi saya tidak tahu bagaimana mengatur sisa parameter dalam fungsi. Bisakah seseorang membantu saya dengan ini?