Cara mensimulasikan analisis kekuatan khusus model lm (menggunakan R)

Mengikuti pertanyaan terakhir yang kami miliki di sini .

Saya ingin tahu apakah ada yang menemukan atau dapat berbagi kode R untuk melakukan analisis kekuatan khusus berdasarkan simulasi untuk model linier?

Kemudian saya jelas ingin memperluasnya ke model yang lebih kompleks, tetapi saya tampaknya tempat yang tepat untuk memulai. Terima kasih.

Saya tidak yakin Anda perlu simulasi untuk model regresi sederhana. Sebagai contoh, lihat kertas Portable Power . Untuk model yang lebih kompleks, khususnya efek campuran, paket pamm di R melakukan analisis daya melalui simulasi. Juga lihat posting Todd Jobe yang memiliki kode R untuk simulasi.

Berikut adalah beberapa sumber kode simulasi dalam R. Saya tidak yakin apakah ada yang secara khusus membahas model linear, tetapi mungkin mereka memberikan cukup contoh untuk mendapatkan intinya:

• Benjamin Bolker telah menulis sebuah buku besar data Ekologis dan Model dengan R . Draf awal seluruh buku bersama dengan kode Sweave tersedia online. Bab 5 membahas analisis daya dan simulasi.

Ada beberapa contoh simulasi di situs berikut:

• http://www.personality-project.org/R/r.datageneration.html
• http://psy-ed.wikidot.com/simulation

Diadaptasi dari Bolker 2009 Model dan Data Ekologi dalam R. Anda perlu menyatakan kekuatan tren (yaitu kemiringan) yang ingin Anda uji. Secara intuitif tren yang kuat dan variabilitas yang rendah akan membutuhkan ukuran sampel yang kecil, tren yang lemah dan variabilitas yang besar akan membutuhkan ukuran sampel yang besar.

a = 2  #desired slope
b = 1  #estimated intercept
sd = 20  #estimated variability defined by standard deviation
nsim = 400  #400 simulations
pval = numeric(nsim)  #placeholder for the second for loop output
Nvec = seq(25, 100, by = 1)  #vector for the range of sample sizes to be tested
power.N = numeric(length(Nvec))   #create placeholder for first for loop output
for (j in 1:length(Nvec)) {
  N = Nvec[j]  
  x = seq(1, 20, length = Nvec[j])  #x value length needs to match sample size (Nvec) length
  for (i in 1:nsim) {   #for this value of N, create random error 400 times
    y_det = a + b * x
    y = rnorm(N, mean = y_det, sd = sd)
    m = lm(y ~ x)
    pval[i] = coef(summary(m))["x", "Pr(>|t|)"]  #all the p values for 400 sims
  }  #cycle through all N values
  power.N[j] = sum(pval < 0.05)/nsim  #the proportion of correct p-values (i.e the power)
}
power.N
plot(Nvec, power.N)  #need about 90 - 100 samples for 80% power

Anda juga dapat mensimulasikan apa tren minimum yang dapat Anda uji untuk ukuran sampel tertentu, seperti yang ditunjukkan dalam buku ini

bvec = seq(-2, 2, by = 0.1)
power.b = numeric(length(bvec))
for (j in 1:length(bvec)) {
  b = bvec[j]
   for (i in 1:nsim) {
     y_det = a + b * x
     y = rnorm(N, mean = y_det, sd = sd)
     m = lm(y ~ x)
     pval[i] = coef(summary(m))["x", "Pr(>|t|)"]
     }
   power.b[j] = sum(pval < 0.05)/nsim
  }
 power.b
 plot(bvec, power.b)