Power analysis of HLM LMM
1.0.0
Karena reproduktifitas psikologi semakin dihargai, permintaan untuk perhitungan daya untuk analisis statistik juga meningkat. Perhitungan daya dari sebagian besar metode statistik dasar ( t test , ANOVA , dll.) Dapat diselesaikan dalam beberapa perangkat lunak (seperti jamovi ). Namun, sebagai metode statistik canggih, perhitungan daya model campuran linier (LMM) tidak dapat dilakukan dalam perangkat lunak seperti jamovi dengan antarmuka visual, dan paket dan fungsi data terkait di R masih diperlukan. Artikel ini terutama memperkenalkan cara melakukan perhitungan daya LMM di R , termasuk: Pengantar data dan fungsinya, skenario yang berlaku dan terbatas, dan contoh demonstrasi spesifik.
Paket data simr adalah paket yang lebih nyaman untuk analisis daya LMM dalam R Dalam analisis, data model disimulasikan oleh simulasi Monte Carlo dan proporsi yang signifikan dihitung dalam jumlah simulasi yang ditetapkan. Fungsi utama adalah powerSim() .
Hitung kekuatan efek posterior, yaitu, menghitung daya berdasarkan efek dalam model yang dibangun berdasarkan data yang ada;
Hitung kekuatan efek apriori, yaitu menghitung daya berdasarkan ukuran efek dalam penelitian yang ada;
Perkirakan jumlah subjek/item, karena ukuran daya dipengaruhi oleh jumlah subjek atau item. Model spesifik adalah bahwa semakin banyak angka, semakin besar daya, yang memberi kita fungsi memperkirakan jumlah subjek/item yang diperlukan untuk percobaan formal berdasarkan efek pra-eksperimen;
Di sini kami menggunakan data dari percobaan psikologi (alamat unduhan data) untuk menunjukkan. Eksperimen ini mengumpulkan reaksi 10 subjek ( subj ) saat melakukan tugas kognitif ( DV ). Eksperimen memiliki dua kondisi ( CondA dan CondB ), yang dirancang untuk 2*2 subjek dan dalam proyek ( item ). Artinya, dua faktor acak dipertimbangkan pada saat yang sama.
Karena kami fokus pada efek tetap dan efek non-acak di sini, kami hanya mempertimbangkan intersep acak pada bagian efek acak ketika pemodelan, dan bagian efek tetap meneliti dua efek utama dan interaksinya.
Pemodelan adalah sebagai berikut:
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item))
Bagian efek tetap dari model adalah:
> summary(Model)$coef
# Estimate Std. Error df t value Pr(>|t|)
# (Intercept) 279.43090 23.37537 11.71977 11.9540740 6.422210e-08
# CondAA2 24.29565 13.49694 498.43142 1.8000866 7.245150e-02
# CondBB2 12.18484 13.36445 509.31328 0.9117351 3.623395e-01
# CondAA2:CondBB2 -32.82881 19.29629 509.76684 -1.7013020 8.949610e-02
Di sini kita mengambil efek tetap dari memeriksa interaksi antara keduanya ( CondAA2:CondBB2 ) sebagai contoh:
> PowerAB_ttest = simr::powerSim(fit = Model, # 要考察的模型
test = fixed('CondAA2:CondBB2', # 要考察的固定效应的名称
method = 't'), # 选取检验方法，因为固定效应为t检验，因此method设置为t
nsim=50) # 设置模拟次数，建议设置为500 (此时可以获取到较稳定的power)
> PowerAB_ttest
# Power for predictor 'CondAA2:CondBB2', (95% confidence interval):
# 46.00% (31.81, 60.68)
#
# Test: t-test with Satterthwaite degrees of freedom (package lmerTest)
# Effect size for CondAA2:CondBB2 is -33.
#
# Based on 50 simulations, (2 warnings, 0 errors)
# alpha = 0.05, nrow = 553
#
# Time elapsed: 0 h 0 m 8 s
#
# nb: result might be an observed power calculation
Kekuatan efek tetap dari interaksi adalah 0.46 (interval kepercayaan 0.318 – 0.607 ).
Demikian pula, daya yang menghitung efek tetap dari CondA dapat diimplementasikan melalui kode berikut:
> PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
Efek utama/bagian interaksi dari model adalah:
> anova(Model)
# Type III Analysis of Variance Table with Satterthwaite's method
# Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
# CondA 8546 8546 1 489.18 0.6694 0.4137
# CondB 2453 2453 1 509.78 0.1921 0.6613
# CondA:CondB 36951 36951 1 509.77 2.8944 0.0895 .
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Di sini kami memeriksa efek utama CondA sebagai contoh:
> PowerA_Ftest = simr::powerSim(fit = Model, # 要考察的模型
test = fixed('CondA', # 要考察的主效应/交互作用的名称
method = 'f'), # 选取检验方法，因为主效应为f检验，因此method设置为f
nsim=50) # 设置模拟次数，建议设置为500 (此时可以获取到较稳定的power)
> PowerA_Ftest
# Power for predictor 'CondA', (95% confidence interval):
# 12.00% ( 4.53, 24.31)
#
# Test: Type-II F-test (package car)
#
# Based on 50 simulations, (50 warnings, 0 errors)
# alpha = 0.05, nrow = 553
#
# Time elapsed: 0 h 0 m 58 s
#
# nb: result might be an observed power calculation
Kekuatan efek utama CondA adalah 0.12 (interval kepercayaan adalah 0.045 – 0.243 ).
Jika ukuran efek yang relatif stabil telah diperoleh dalam penelitian sebelumnya, misalnya, ditemukan bahwa efek tetap CondA umumnya sekitar 50, daya juga dapat dihitung berdasarkan efek sebelumnya. Langkah -langkahnya adalah sebagai berikut:
> fixef(Model) # 查看模型固定效应
# (Intercept) CondAA2 CondBB2 CondAA2:CondBB2
# 279.43090 24.29565 12.18484 -32.82881
> fixef(Model)[2] = 50 # 修改CondA的效应
> fixef(Model) # 查看修改后的固定效应
# (Intercept) CondAA2 CondBB2 CondAA2:CondBB2
# 279.43090 50.00000 12.18484 -32.82881
> PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
> PowerA_ttest
# Power for predictor 'CondAA2', (95% confidence interval):
# 88.00% (75.69, 95.47)
#
# Test: t-test with Satterthwaite degrees of freedom (package lmerTest)
# Effect size for CondAA2 is 50.
#
# Based on 50 simulations, (3 warnings, 0 errors)
# alpha = 0.05, nrow = 553
#
# Time elapsed: 0 h 0 m 9 s
Pada saat ini, kekuatan efek tetap CondA adalah 0.88 .
> PowerA_Ftest = simr::powerSim(fit = Model, test = fixed('CondA', method = 'f'), nsim=50)
> PowerA_Ftest
# Power for predictor 'CondA', (95% confidence interval):
# 94.00% (83.45, 98.75)
#
# Test: Type-II F-test (package car)
#
# Based on 50 simulations, (50 warnings, 0 errors)
# alpha = 0.05, nrow = 553
#
# Time elapsed: 0 h 1 m 4 s
Efek utama kekuatan CondA adalah 0.94 .
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item)) # 这里仍以后验power为例
Kekuatan efek tetap dari interaksi adalah 0.46 , dan kami ingin menguji jumlah subjek, daya dapat mencapai 0.8 . Langkah -langkahnya adalah sebagai berikut:
> Model2 = extend(object = Model, along = 'subj', n = 30)
> PowerA_ttest2 = powerSim(Model2, fixed('CondAA2', method = 't'), nsim=50)
> PowerA_ttest2
# Power for predictor 'CondAA2', (95% confidence interval):
# 84.00% (70.89, 92.83)
#
# Test: t-test with Satterthwaite degrees of freedom (package lmerTest)
# Effect size for CondAA2 is 24.
#
# Based on 50 simulations, (6 warnings, 0 errors)
# alpha = 0.05, nrow = 1659
#
# Time elapsed: 0 h 0 m 11 s
#
# nb: result might be an observed power calculation
Pada saat ini, daya telah naik menjadi 0.84 .
powerCurve() > Pcurve = powerCurve(fit = Model2, test = fixed('CondAA2', method = 't'), along = 'subj', nsim=50)
> Pcurve
# Power for predictor 'CondAA2', (95% confidence interval), by largest value of subj:
# 11: 22.00% (11.53, 35.96) - 180 rows
# 14: 36.00% (22.92, 50.81) - 325 rows
# 17: 50.00% (35.53, 64.47) - 507 rows
# 2: 58.00% (43.21, 71.81) - 667 rows
# 22: 64.00% (49.19, 77.08) - 837 rows
# 25: 58.00% (43.21, 71.81) - 989 rows
# 28: 74.00% (59.66, 85.37) - 1166 rows
# 30: 82.00% (68.56, 91.42) - 1318 rows
# 6: 86.00% (73.26, 94.18) - 1489 rows
# 9: 90.00% (78.19, 96.67) - 1659 rows
#
# Time elapsed: 0 h 1 m 39 s
Perhatikan bahwa nilai di depan setiap daya tidak akurat saat ini, dan nilai sebenarnya adalah:
> Pcurve$nlevels
# [1] 3 6 9 12 15 18 21 24 27 30
Dapat dilihat bahwa ketika sekitar 24 subjek mencapai daya, daya dapat mencapai lebih dari 0.8 .
Kisaran ukuran sampel dapat diatur secara manual:
> Pcurve = powerCurve(fit = Model2, test = fixed('CondAA2', method = 't'), along = 'subj', nsim=50,
breaks=c(21, 22, 23, 24, 25, 30)) # 考察样本量为21、22、23、24、25、30时的power
> Pcurve
# Power for predictor 'CondAA2', (95% confidence interval), by largest value of subj:
# 28: 72.00% (57.51, 83.77) - 1166 rows
# 29: 74.00% (59.66, 85.37) - 1220 rows
# 3: 70.00% (55.39, 82.14) - 1262 rows
# 30: 74.00% (59.66, 85.37) - 1318 rows
# 4: 76.00% (61.83, 86.94) - 1369 rows
# 9: 80.00% (66.28, 89.97) - 1659 rows
#
# Time elapsed: 0 h 1 m 12 s
> Pcurve$nlevels
# [1] 21 22 23 24 25 30
Perhitungan daya efek dalam model umum mirip dengan yang di atas, perbedaan utamanya adalah:
method dalam fixed() harus diatur ke 'z' karena efek tetap dalam model umum adalah uji-z;method dalam fixed() harus diatur ke 'chisq' , karena efek utama dalam model umum adalah uji chi-square
