Power analysis of HLM LMM
1.0.0
Comme la reproductibilité de la psychologie est de plus en plus valorisée, la demande de calculs de pouvoir pour l'analyse statistique augmente également. Les calculs de puissance de la plupart des méthodes statistiques de base ( t test , ANOVA , etc.) peuvent être achevés dans certains logiciels (comme jamovi ). Cependant, en tant que méthode statistique avancée, le calcul de puissance du modèle mixte linéaire (LMM) ne peut pas être effectué dans des logiciels tels que jamovi avec interface visuelle, et des paquets de données et des fonctions connexes en R sont toujours nécessaires. Cet article présente principalement comment effectuer le calcul de la puissance de LMM dans R , notamment: Introduction aux données et ses fonctions, scénarios applicables et restreints, et des exemples de démonstrations spécifiques.
simr est un package plus pratique pour l'analyse de la puissance LMM dans R Dans l'analyse, les données du modèle sont simulées par la simulation Monte Carlo et la proportion significative est calculée dans le nombre d'ensembles de simulations. La fonction principale est powerSim() .
Calculez la puissance de l'effet postérieur, c'est-à-dire calculer la puissance en fonction des effets du modèle construits en fonction des données existantes;
Calculez la puissance de l'effet a priori, c'est-à-dire calculer la puissance en fonction de la taille d'un effet dans la recherche existante;
Estimez le nombre de sujets / éléments, car la taille de la puissance est affectée par le nombre de sujets ou d'éléments. Le modèle spécifique est que plus il y a de nombre, plus la puissance est grande, ce qui nous offre la fonction d'estimation du nombre de sujets / éléments requis pour les expériences formelles en fonction des effets de la pré-expérience;
Ici, nous utilisons les données d'une expérience de psychologie (adresse de téléchargement des données) pour démontrer. Cette expérience a collecté les réactions de 10 sujets ( subj ) lors de l'exécution d'une tâche cognitive ( DV ). L'expérience avait deux conditions ( CondA et CondB ), qui ont été conçues pour 2*2 et dans le projet ( item ). Autrement dit, deux facteurs aléatoires sont pris en compte en même temps.
Parce que nous nous concentrons sur les effets fixes et les effets non alétiques ici, nous ne considérons que des interceptions aléatoires dans la partie des effets aléatoires lors de la modélisation, et la partie des effets fixes examine les deux effets principaux et leurs interactions.
La modélisation est la suivante:
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item))
La partie à effet fixe du modèle est:
> 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
Ici, nous prenons l'effet fixe de l'examen de l'interaction entre les deux ( CondAA2:CondBB2 ) à titre d'exemple:
> 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
La puissance de l'effet fixe de l'interaction était 0.46 (intervalle de confiance 0.318 – 0.607 ).
De même, la puissance qui calcule l'effet fixe de CondA peut être implémentée via le code suivant:
> PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
La partie principale de l'effet / interaction du modèle est:
> 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
Ici, nous prenons en examinant l'effet principal de CondA à titre d'exemple:
> 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
La puissance de l'effet principal de CondA est 0.12 (l'intervalle de confiance est 0.045 – 0.243 ).
Si une taille d'effet relativement stable a été obtenue dans des études précédentes, par exemple, il est constaté que l'effet fixe de CondA est généralement d'environ 50, la puissance peut également être calculée en fonction des effets antérieurs. Les étapes sont les suivantes:
> 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
À l'heure actuelle, la puissance à effet fixe de CondA est 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
La puissance d'effet principale de CondA est 0.94 .
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item)) # 这里仍以后验power为例
La puissance de l'effet fixe de l'interaction est 0.46 et nous voulons tester le nombre de sujets, la puissance peut atteindre 0.8 . Les étapes sont les suivantes:
> 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
À l'heure actuelle, la puissance est passée à 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
Notez que la valeur devant chaque puissance est inexacte pour le moment, et la valeur réelle est:
> Pcurve$nlevels
# [1] 3 6 9 12 15 18 21 24 27 30
On peut voir qu'environ 24 sujets atteignent la puissance, la puissance peut atteindre plus de 0.8 .
La plage de taille de l'échantillon peut être réglée manuellement:
> 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
Le calcul de puissance de l'effet dans le modèle généralisé est similaire à ce qui précède, les principales différences sont:
method dans fixed() doit être défini sur 'z' car l'effet fixe dans le modèle généralisé est le test Z;method dans fixed() doit être défini sur 'chisq' , car l'effet principal dans le modèle généralisé est un test du chi carré
