Power analysis of HLM LMM
1.0.0
เนื่องจากความสามารถในการทำซ้ำของจิตวิทยามีค่ามากขึ้นเรื่อย ๆ ความต้องการการคำนวณพลังงานสำหรับการวิเคราะห์ทางสถิติก็เพิ่มขึ้นเช่นกัน การคำนวณพลังงานของวิธีการทางสถิติพื้นฐานที่สุด ( t test , ANOVA ฯลฯ ) สามารถทำได้ในซอฟต์แวร์บางอย่าง (เช่น jamovi ) อย่างไรก็ตามในฐานะที่เป็นวิธีการทางสถิติขั้นสูงการคำนวณพลังงานของแบบจำลองเชิงเส้นแบบผสม (LMM) ไม่สามารถทำได้ในซอฟต์แวร์เช่น jamovi กับอินเทอร์เฟซภาพและแพ็คเก็ตข้อมูลและฟังก์ชั่นที่เกี่ยวข้องใน R ยังคงเป็นสิ่งจำเป็น บทความนี้ส่วนใหญ่แนะนำวิธีการคำนวณพลังงานของ LMM ใน R รวมถึง: บทนำสู่ข้อมูลและฟังก์ชั่นของมันสถานการณ์ที่ใช้งานได้และ จำกัด และการสาธิตตัวอย่างที่เฉพาะเจาะจง
แพ็คเก็ต Data simr เป็นแพ็คเกจที่สะดวกกว่าสำหรับการวิเคราะห์พลังงาน LMM ใน R ในการวิเคราะห์ข้อมูลแบบจำลองจะถูกจำลองโดยการจำลอง Monte Carlo และสัดส่วนที่สำคัญคำนวณในจำนวนชุดของการจำลอง ฟังก์ชั่นหลักคือ powerSim()
คำนวณกำลังของเอฟเฟกต์หลังนั่นคือคำนวณพลังงานตามเอฟเฟกต์ในแบบจำลองที่สร้างขึ้นตามข้อมูลที่มีอยู่
คำนวณพลังของเอฟเฟกต์เบื้องต้นนั่นคือคำนวณพลังงานตามขนาดของผลกระทบในการวิจัยที่มีอยู่
ประเมินจำนวนวิชา/รายการเนื่องจากขนาดของพลังงานได้รับผลกระทบจากจำนวนวิชาหรือรายการ โมเดลเฉพาะคือยิ่งมีจำนวนมากเท่าใดก็ยิ่งมีขนาดใหญ่ขึ้นซึ่งจะช่วยให้เราได้ฟังก์ชั่นการประมาณจำนวนของวิชา/รายการที่จำเป็นสำหรับการทดลองอย่างเป็นทางการตามผลกระทบของการทดลองก่อน
ที่นี่เราใช้ข้อมูลจากการทดลองทางจิตวิทยา (ที่อยู่ดาวน์โหลดข้อมูล) เพื่อสาธิต การทดลองนี้รวบรวมปฏิกิริยาของอาสาสมัคร 10 คน ( subj ) เมื่อดำเนินงานเกี่ยวกับความรู้ความเข้าใจ ( DV ) การทดลองมีสองเงื่อนไข ( CondA และ CondB ) ซึ่งออกแบบมาสำหรับ 2*2 วิชาและภายในโครงการ ( item ) นั่นคือสองปัจจัยสุ่มได้รับการพิจารณาในเวลาเดียวกัน
เนื่องจากเรามุ่งเน้นไปที่เอฟเฟกต์คงที่และเอฟเฟกต์ที่ไม่สุ่มที่นี่เราจึงพิจารณาการสกัดกั้นแบบสุ่มในส่วนเอฟเฟกต์แบบสุ่มเมื่อการสร้างแบบจำลองและส่วนเอฟเฟกต์คงที่จะตรวจสอบผลกระทบหลักทั้งสองและการโต้ตอบของพวกเขา
การสร้างแบบจำลองมีดังนี้:
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item))
ส่วนเอฟเฟกต์คงที่ของโมเดลคือ:
> 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
ที่นี่เรามีผลคงที่ของการตรวจสอบปฏิสัมพันธ์ระหว่างทั้งสอง ( CondAA2:CondBB2 ) เป็นตัวอย่าง:
> 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
พลังของผลคงที่ของการโต้ตอบคือ 0.46 (ช่วงความเชื่อมั่น 0.318 – 0.607 )
ในทำนองเดียวกันพลังงานที่คำนวณผลคงที่ของ CondA สามารถนำไปใช้ผ่านรหัสต่อไปนี้:
> PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
ผลกระทบหลัก/ส่วนการโต้ตอบของโมเดลคือ:
> 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
ที่นี่เราตรวจสอบผลกระทบหลักของ CondA เป็นตัวอย่าง:
> 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
พลังของผลกระทบหลักของ CondA คือ 0.12 (ช่วงความเชื่อมั่นคือ 0.045 – 0.243 )
หากขนาดของเอฟเฟกต์ที่ค่อนข้างเสถียรได้รับในการศึกษาก่อนหน้านี้ตัวอย่างเช่นพบว่าผลคงที่ของ CondA โดยทั่วไปประมาณ 50 พลังงานสามารถคำนวณได้ตามผลกระทบก่อนหน้านี้ ขั้นตอนมีดังนี้:
> 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
ในเวลานี้พลังเอฟเฟกต์คงที่ของ CondA คือ 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
พลังหลักของ CondA คือ 0.94
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item)) # 这里仍以后验power为例
พลังของเอฟเฟกต์คงที่ของการโต้ตอบคือ 0.46 และเราต้องการทดสอบจำนวนวิชาพลังสามารถถึง 0.8 ขั้นตอนมีดังนี้:
> 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
ในเวลานี้พลังงานเพิ่มขึ้นเป็น 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
โปรดทราบ ว่าค่าที่อยู่ด้านหน้าของพลังงานแต่ละอันไม่ถูกต้องในเวลานี้และค่าที่แท้จริงคือ:
> Pcurve$nlevels
# [1] 3 6 9 12 15 18 21 24 27 30
จะเห็นได้ว่าเมื่อประมาณ 24 วิชาถึงพลังอำนาจสามารถเข้าถึงได้มากกว่า 0.8
ช่วงขนาดตัวอย่างสามารถตั้งค่าด้วยตนเอง:
> 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
การคำนวณพลังงานของเอฟเฟกต์ในแบบจำลองทั่วไปนั้นคล้ายกับข้างต้นความแตกต่างหลักคือ:
method ใน fixed() ควรตั้งค่าเป็น 'z' เนื่องจากเอฟเฟกต์คงที่ในโมเดลทั่วไปคือการทดสอบ Z;method ใน fixed() ควรตั้งค่าเป็น 'chisq' เนื่องจากเอฟเฟกต์หลักในแบบจำลองทั่วไปคือการทดสอบไคสแควร์
