Power analysis of HLM LMM
1.0.0
Поскольку воспроизводимость психологии все больше ценится, спрос на расчеты власти для статистического анализа также увеличивается. Расчеты мощности самых основных статистических методов ( t test , ANOVA и т. Д.) Могут быть завершены в некотором программном обеспечении (например, jamovi ). Однако в качестве расширенного статистического метода расчет мощности линейной смешанной модели (LMM) не может быть выполнен в программном обеспечении, таком как jamovi , с визуальным интерфейсом, а соответствующие пакеты данных и функции в R в R все еще требуются. Эта статья в основном представляет, как выполнить расчет мощности LMM в R , в том числе: Введение в данные и его функции, применимые и ограниченные сценарии и конкретные примеры демонстрации.
Пакет данных simr - это пакет, который более удобен для анализа мощности LMM в R В анализе данные модели моделируются моделированием Монте -Карло, и значительная доля рассчитывается в определенном количестве моделирования. Основной функцией является 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' , потому что основным эффектом в обобщенной модели является тест хи-квадрат
