Pertanyaan ini sebagai jawaban atas jawaban yang diberikan oleh @Greg Snow sehubungan dengan pertanyaan yang saya ajukan mengenai analisis daya dengan regresi logistik dan SAS Proc GLMPOWER.

Jika saya merancang percobaan dan akan menganalisis hasil dalam regresi logistik faktorial, bagaimana saya bisa menggunakan simulasi (dan di sini ) untuk melakukan analisis daya?

Berikut adalah contoh sederhana di mana ada dua variabel, yang pertama mengambil tiga nilai yang mungkin {0,03, 0,06, 0,09} dan yang kedua adalah indikator dummy {0,1}. Untuk masing-masing kami perkirakan tingkat respons untuk setiap kombinasi (# responden / jumlah orang yang dipasarkan). Lebih lanjut, kami ingin memiliki 3 kali lebih banyak kombinasi faktor pertama dari yang lainnya (yang dapat dianggap sama) karena kombinasi pertama ini adalah versi kami yang telah dicoba dan benar. Ini adalah pengaturan seperti yang diberikan dalam kursus SAS yang disebutkan dalam pertanyaan terkait.

Model yang akan digunakan untuk menganalisis hasil adalah regresi logistik, dengan efek dan interaksi utama (responsnya 0 atau 1).

mod <- glm(response ~ Var1 + Var2 + I(Var1*Var2))

Bagaimana saya bisa mensimulasikan set data untuk digunakan dengan model ini untuk melakukan analisis daya?

Ketika saya menjalankan ini melalui SAS Proc GLMPOWER(menggunakan STDDEV =0.05486016 yang sesuai dengan di sqrt(p(1-p))mana p adalah rata-rata tertimbang dari tingkat respons yang ditampilkan):

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;

Catatan: GLMPOWERhanya akan menggunakan variabel kelas (nominal) sehingga 3, 6, 9 di atas diperlakukan sebagai karakter dan bisa saja rendah, sedang dan tinggi atau tiga string lainnya. Ketika analisis sebenarnya dilakukan, Var1 akan digunakan numerik (dan kami akan menyertakan istilah polinomial Var1 * Var1) untuk menjelaskan setiap kelengkungan.

Output dari SAS adalah

Jadi kita melihat bahwa kita perlu 762.112 sebagai ukuran sampel kami (efek utama Var2 adalah yang paling sulit untuk diperkirakan) dengan kekuatan sama dengan 0,80 dan alpha sama dengan 0,05. Kami akan mengalokasikan ini sehingga 3 kali lebih banyak dari kombinasi baseline (yaitu 0,375 * 762112) dan sisanya hanya jatuh secara merata ke dalam 5 kombinasi lainnya.

Persiapan:

• Seperti dibahas dalam manual G * Power , ada beberapa jenis analisis daya, tergantung pada apa yang ingin Anda selesaikan. (Yaitu, , ukuran efek , , dan power ada dalam hubungannya satu sama lain; dengan menentukan tiga dari mereka akan membiarkan Anda menyelesaikan untuk yang keempat.) E S $N$$ES$$\alpha$

  • dalam deskripsi Anda, Anda ingin tahu sesuai untuk menangkap tingkat respons yang Anda tentukan dengan , dan daya = 80%. Ini adalah kekuatan a-priori . α = $N$$\alpha=.05$
  • kita dapat mulai dengan kekuatan post-hoc (tentukan daya yang diberikan , tingkat respons, & alpha) karena ini secara konsep lebih sederhana, dan kemudian naik$N$

• Selain postingan @ GregSnow yang luar biasa , panduan hebat lainnya untuk analisis daya berbasis simulasi pada CV dapat ditemukan di sini: Menghitung kekuatan statistik . Untuk merangkum ide-ide dasar:

  1. mencari tahu efek yang Anda ingin dapat dideteksi
  2. menghasilkan data N dari dunia yang mungkin
  3. jalankan analisis yang ingin Anda lakukan atas data palsu itu
  4. simpan apakah hasilnya 'signifikan' sesuai dengan alpha yang Anda pilih
  5. ulangi banyak ( ) kali & gunakan% 'signifikan' sebagai perkiraan daya (post-hoc) di$B$$N$
  6. untuk menentukan kekuatan a-priori, cari kemungkinan untuk menemukan nilai yang menghasilkan kekuatan yang Anda inginkan $N$

• Apakah Anda akan menemukan signifikansi pada iterasi tertentu dapat dipahami sebagai hasil dari percobaan Bernoulli dengan probabilitas (di mana adalah kekuatan). Proporsi yang ditemukan di atas iterasi memungkinkan kita untuk memperkirakan sebenarnya . Untuk mendapatkan perkiraan yang lebih baik, kita dapat meningkatkan , meskipun ini juga akan membuat simulasi lebih lama. p B p $p$$p$$B$$p$$B$

• Dalam R, cara utama untuk menghasilkan data biner dengan probabilitas keberhasilan yang diberikan adalah ? Rbinom

  • Misalnya untuk mendapatkan jumlah keberhasilan dari 10 percobaan Bernoulli dengan probabilitas p, kodenya adalah rbinom(n=10, size=1, prob=p), (Anda mungkin ingin menugaskan hasilnya ke variabel untuk penyimpanan)
  • Anda juga dapat menghasilkan data seperti itu dengan kurang elegan dengan menggunakan ? runif , misalnya,ifelse(runif(1)<=p, 1, 0)
  • jika Anda yakin hasilnya dimediasi oleh variabel Gaussian laten, Anda dapat menghasilkan variabel laten sebagai fungsi kovariat Anda dengan ? rnorm , dan kemudian mengubahnya menjadi probabilitas dengan pnorm()dan menggunakannya dalam rbinom()kode Anda .

• Anda menyatakan bahwa Anda akan "memasukkan istilah polinomial Var1 * Var1) untuk memperhitungkan setiap kelengkungan". Ada kebingungan di sini; istilah polinom dapat membantu kita menjelaskan kelengkungan, tetapi ini adalah istilah interaksi - itu tidak akan membantu kita dengan cara ini. Meskipun demikian, tingkat respons Anda mengharuskan kami untuk memasukkan istilah kuadrat dan istilah interaksi dalam model kami. Khususnya, model Anda perlu menyertakan: , , dan , di luar ketentuan dasar. v a r 1 ∗ v a r 2 v a r 1 2 ∗ v a r $var1^2$$var1*var2$$var1^2*var2$

• Meskipun ditulis dalam konteks pertanyaan yang berbeda, jawaban saya di sini: Perbedaan antara model logit dan probit memiliki banyak informasi dasar tentang jenis model ini.

• Seperti halnya ada berbagai jenis Tipe I tingkat kesalahan ketika ada beberapa hipotesis (misalnya, tingkat kesalahan per-kontras , tingkat kesalahan familywise , & tingkat kesalahan per-keluarga ), sehingga ada berbagai jenis kekuasaan (misalnya, untuk sebuah efek yang ditentukan sebelumnya , untuk efek apa pun , & untuk semua efek ). Anda juga dapat mencari kekuatan untuk mendeteksi kombinasi efek tertentu, atau untuk kekuatan uji simultan model secara keseluruhan. Dugaan saya dari deskripsi kode SAS Anda adalah bahwa ia mencari yang terakhir. Namun, dari uraian situasi Anda, saya mengasumsikan Anda ingin mendeteksi efek interaksi seminimal mungkin.

  • * referensi: Maxwell, SE (2004). Kegigihan studi yang kurang kuat dalam penelitian psikologis: penyebab, konsekuensi, dan perbaikan. Metode Psikologis , 9 , 2 , hlm. 147-163.
  • efek Anda cukup kecil (jangan dikacaukan dengan tingkat respons yang rendah), jadi kami akan kesulitan mencapai kekuatan yang baik.
  • Perhatikan bahwa, meskipun semua ini terdengar sangat mirip, mereka sangat tidak sama (misalnya, sangat mungkin untuk mendapatkan model yang signifikan tanpa efek yang signifikan - dibahas di sini: Bagaimana regresi menjadi signifikan namun semua prediktor tidak signifikan?, atau efek signifikan tetapi di mana model tidak signifikan - dibahas di sini: Signifikansi koefisien dalam regresi linier: t-test signifikan vs F-statistik tidak signifikan ), yang akan diilustrasikan di bawah ini.

• Untuk cara berbeda dalam memikirkan masalah yang berkaitan dengan kekuasaan, lihat jawaban saya di sini: Cara melaporkan ketepatan umum dalam memperkirakan korelasi dalam konteks justifikasi ukuran sampel.

Kekuatan post-hoc sederhana untuk regresi logistik di R:

Katakanlah tingkat respons yang Anda ajukan mewakili situasi sebenarnya di dunia, dan Anda telah mengirim 10.000 surat. Apa kekuatan untuk mendeteksi efek-efek itu? (Perhatikan bahwa saya terkenal karena menulis kode "komikal tidak efisien", yang berikut ini dimaksudkan untuk mudah diikuti daripada dioptimalkan untuk efisiensi; pada kenyataannya, ini sangat lambat.)

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

Jadi kita melihat bahwa 10.000 huruf tidak benar-benar mencapai daya 80% (dalam bentuk apa pun) untuk mendeteksi tingkat respons ini. (Saya tidak cukup yakin tentang apa yang dilakukan kode SAS untuk dapat menjelaskan perbedaan mencolok antara pendekatan ini, tetapi kode ini secara konseptual mudah - jika lambat - dan saya telah meluangkan waktu untuk memeriksanya, dan saya pikir ini hasilnya masuk akal.)

Kekuatan a-priori berbasis simulasi untuk regresi logistik:

Dari sini idenya adalah hanya untuk mencari kemungkinan sampai kita menemukan nilai yang menghasilkan tingkat yang diinginkan dari jenis daya yang Anda minati. Setiap strategi pencarian yang dapat Anda kodekan untuk bekerja dengan ini akan baik-baik saja (dalam teori). Mengingat 's yang akan diperlukan untuk menangkap efek kecil seperti itu, ada baiknya memikirkan bagaimana melakukan ini dengan lebih efisien. Pendekatan khas saya hanyalah kekuatan kasar, yaitu untuk menilai setiap yang mungkin saya pertimbangkan. (Namun perhatikan, bahwa saya biasanya hanya mempertimbangkan kisaran kecil, dan saya biasanya bekerja dengan sangat kecil - setidaknya dibandingkan dengan ini.) N N $N$$N$$N$$N$

Alih-alih, strategi saya di sini adalah mengurung kemungkinan untuk mengetahui kisaran daya yang akan terjadi. Jadi, saya mengambil dari 500.000 dan menjalankan kembali kode (memulai benih yang sama, nb ini membutuhkan waktu satu setengah jam untuk menjalankan). Inilah hasilnya: $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

Kita dapat melihat dari sini bahwa besarnya efek Anda sangat bervariasi, dan dengan demikian kemampuan Anda untuk mendeteksinya bervariasi. Misalnya, efek sangat sulit dideteksi, hanya signifikan 6% dari waktu bahkan dengan setengah juta huruf. Di sisi lain, model secara keseluruhan selalu jauh lebih baik daripada model nol. Kemungkinan lain tersusun di antara keduanya. Meskipun sebagian besar 'data' dibuang pada setiap iterasi, sedikit eksplorasi masih mungkin dilakukan. Sebagai contoh, kita dapat menggunakan matriks untuk menilai korelasi antara probabilitas variabel yang berbeda menjadi signifikan. $var1^2$significant

Saya harus mencatat sebagai kesimpulan, bahwa karena kompleksitas dan besar yang diperlukan dalam situasi Anda, ini tidak sesederhana seperti yang saya duga / klaim dalam komentar awal saya. Namun, Anda tentu bisa mendapatkan ide untuk bagaimana hal ini dapat dilakukan secara umum, dan masalah yang terlibat dalam analisis kekuatan, dari apa yang saya taruh di sini. HTH. $N$

Jawaban @ Gung sangat bagus untuk pengertian. Inilah pendekatan yang akan saya gunakan:

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

Fungsi ini menguji efek keseluruhan dari v2, model dapat diubah untuk melihat jenis tes lainnya. Saya suka menulisnya sebagai fungsi sehingga ketika saya ingin menguji sesuatu yang berbeda saya hanya bisa mengubah argumen fungsi. Anda juga bisa menambahkan bilah kemajuan, atau menggunakan paket paralel untuk mempercepat.

Di sini saya hanya melakukan 100 replikasi, saya biasanya mulai sekitar tingkat itu untuk menemukan perkiraan ukuran sampel, kemudian naik iterasi ketika saya berada di taman bola yang tepat (tidak perlu membuang waktu pada 10.000 iterasi ketika Anda memiliki daya 20%).