Power analysis of HLM LMM
1.0.0
Da die Reproduzierbarkeit der Psychologie zunehmend geschätzt wird, nimmt auch die Nachfrage nach Leistungsberechnungen für die statistische Analyse zu. Leistungsberechnungen der meisten grundlegenden statistischen Methoden ( t test , ANOVA usw.) können in einer Software (wie jamovi ) abgeschlossen werden. Als erweiterte statistische Methode kann jedoch die Leistungsberechnung des linearen gemischten Modells (LMM) in Software wie jamovi mit visueller Schnittstelle nicht durchgeführt werden, und verwandte Datenpakete und Funktionen in R sind weiterhin erforderlich. In diesem Artikel wird hauptsächlich eingeführt, wie die Leistungsberechnung von LMM in R durchgeführt wird, einschließlich: Einführung in Daten und ihre Funktionen, anwendbaren und eingeschränkten Szenarien sowie spezifische Beispieldemonstrationen.
simr -Datenpaket ist ein Paket, das für die LMM -Leistungsanalyse in R in der Analyse bequemer ist. In der Analyse werden die Modelldaten durch Monte -Carlo -Simulation simuliert und der erhebliche Anteil wird in der festgelegten Anzahl von Simulationen berechnet. Die Hauptfunktion ist powerSim() .
Berechnen Sie die Leistung des hinteren Effekts, dh die Leistung basierend auf den Auswirkungen des Modells, das auf vorhandenen Daten basiert.
Berechnen Sie die Leistung eines priori -Effekts, dh die Leistung basierend auf der Größe eines Effekts in der vorhandenen Forschung;
Schätzen Sie die Anzahl der Probanden/Elemente, da die Leistungsgröße durch die Anzahl der Probanden oder Elemente beeinflusst wird. Das spezifische Modell ist, dass je mehr Anzahl ist, desto größer ist die Leistung, die uns die Funktion zur Schätzung der Anzahl der Probanden/Elemente bietet, die für formale Experimente erforderlich sind, basierend auf den Auswirkungen des Voraussagens;
Hier verwenden wir die Daten aus einem Psychologie -Experiment (Daten -Download -Adresse), um zu demonstrieren. Dieses Experiment sammelte die Reaktionen von 10 Probanden ( subj ) bei der Ausführung einer kognitiven Aufgabe ( DV ). Das Experiment hatte zwei Bedingungen ( CondA und CondB ), die für 2*2 Probanden und innerhalb des Projekts ( item ) ausgelegt waren. Das heißt, zwei zufällige Faktoren werden gleichzeitig berücksichtigt.
Da wir uns hier auf feste Effekte und nicht zufällige Effekte konzentrieren, betrachten wir nur zufällige Abschnitte im Teil der zufälligen Effekte bei der Modellierung, und der Teil der festen Effekte untersucht die beiden Haupteffekte und deren Wechselwirkungen.
Die Modellierung ist wie folgt:
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item))
Der feste Effektsteil des Modells ist:
> 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
Hier haben wir den festen Effekt, die Wechselwirkung zwischen den beiden ( CondAA2:CondBB2 ) als Beispiel zu untersuchen:
> 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
Die Leistung des festen Effekts der Wechselwirkung betrug 0.46 (Konfidenzintervall 0.318 – 0.607 ).
In ähnlicher Weise kann eine Leistung, die den festen Effekt von CondA berechnet, über den folgenden Code implementiert werden:
> PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
Der Haupteffekt-/Interaktionsteil des Modells ist:
> 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
Hier untersuchen wir den Haupteffekt von CondA als Beispiel:
> 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
Die Leistung des Haupteffekts von CondA beträgt 0.12 (Konfidenzintervall beträgt 0.045 – 0.243 ).
Wenn in früheren Studien eine relativ stabile Effektgröße erhalten wurde, wird beispielsweise festgestellt, dass der feste Effekt von CondA im Allgemeinen etwa 50 beträgt, sondern auch auf der Grundlage der vorherigen Effekte berechnet werden kann. Die Schritte sind wie folgt:
> 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
Zu diesem Zeitpunkt beträgt die feste Effektleistung von 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
Die Haupteffektleistung von CondA beträgt 0.94 .
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item)) # 这里仍以后验power为例
Die Leistung des festen Effekts der Wechselwirkung beträgt 0.46 , und wir möchten die Anzahl der Probanden testen, die Leistung kann 0.8 erreichen. Die Schritte sind wie folgt:
> 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
Zu diesem Zeitpunkt ist die Macht auf 0.84 gestiegen.
powerCurve() verwenden > 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
Beachten Sie , dass der Wert vor jeder Leistung zu diesem Zeitpunkt ungenau ist und der wahre Wert lautet:
> Pcurve$nlevels
# [1] 3 6 9 12 15 18 21 24 27 30
Es ist ersichtlich, dass bei etwa 24 Probanden die Leistung mehr als 0.8 erreichen kann.
Der Stichprobengrößenbereich kann manuell eingestellt werden:
> 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
Die Leistungsberechnung des Effekts im verallgemeinerten Modell ähnelt dem oben. Die Hauptunterschiede sind:
method in fixed() auf 'z' eingestellt werden, da der feste Effekt im verallgemeinerten Modell der Z-Test ist.method in fixed() auf 'chisq' eingestellt werden, da der Haupteffekt im verallgemeinerten Modell ein Chi-Quadrat-Test ist
