Power analysis of HLM LMM
1.0.0
As the reproducibility of psychology is increasingly valued, the demand for power calculations for statistical analysis is also increasing. Power calculations of most basic statistical methods ( t test , ANOVA , etc.) can be completed in some software (such as jamovi ). However, as an advanced statistical method, the power calculation of linear mixed model (LMM) cannot be performed in software such as jamovi with visual interface, and related data packets and functions in R are still required. This article mainly introduces how to perform power calculation of LMM in R , including: introduction to data and its functions, applicable and restricted scenarios, and specific example demonstrations.
simr data packet is a package that is more convenient for LMM power analysis in R In the analysis, the model data is simulated by Monte Carlo simulation and the significant proportion is calculated in the set number of simulations. The main function is powerSim() .
Calculate the power of the posterior effect, that is, calculate the power based on the effects in the model built based on existing data;
Calculate the power of a priori effect, that is, calculate the power based on the size of an effect in existing research;
Estimate the number of subjects/items, because the size of power is affected by the number of subjects or items. The specific model is that the more number, the larger the power is, which provides us with the function of estimating the number of subjects/items required for formal experiments based on the effects of the pre-experiment;
Here we use the data from a psychology experiment (data download address) to demonstrate. This experiment collected the reactions of 10 subjects ( subj ) when performing a cognitive task ( DV ). The experiment had two conditions ( CondA and CondB ), which were designed for 2*2 subjects and within the project ( item ). That is, two random factors are considered at the same time.
Because we focus on fixed effects and non-random effects here, we only consider random intercepts in the random effects part when modeling, and the fixed effects part examines the two main effects and their interactions.
Modeling is as follows:
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item))
The fixed effect part of the model is:
> 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
Here we take the fixed effect of examining the interaction between the two ( CondAA2:CondBB2 ) as an example:
> 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
The power of the fixed effect of the interaction was 0.46 (confidence interval 0.318 – 0.607 ).
Similarly, power that calculates the fixed effect of CondA can be implemented through the following code:
> PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
The main effect/interaction part of the model is:
> 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
Here we take examining the main effect of CondA as an example:
> 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
The power of the main effect of CondA is 0.12 (confidence interval is 0.045 – 0.243 ).
If a relatively stable effect size has been obtained in previous studies, for example, it is found that the fixed effect of CondA is generally around 50, power can also be calculated based on the prior effects. The steps are as follows:
> 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
At this time, the fixed effect power of CondA is 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
The main effect power of CondA is 0.94 .
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item)) # 这里仍以后验power为例
The power of the fixed effect of the interaction is 0.46 , and we want to test the number of subjects, the power can reach 0.8 . The steps are as follows:
> 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
At this time, power has risen to 0.84 .
powerCurve() function > 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
Note that the value in front of each power is inaccurate at this time, and the real value is:
> Pcurve$nlevels
# [1] 3 6 9 12 15 18 21 24 27 30
It can be seen that when about 24 subjects reach the power, the power can reach more than 0.8 .
The sample size range can be set manually:
> 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
The power calculation of the effect in the generalized model is similar to the above, the main differences are:
method parameter in fixed() should be set to 'z' because the fixed effect in the generalized model is the z-test;method parameter in fixed() should be set to 'chisq' , because the main effect in the generalized model is a chi-square test
