Power analysis of HLM LMM
1.0.0
Como a reprodutibilidade da psicologia é cada vez mais valorizada, a demanda por cálculos de energia para análise estatística também está aumentando. Os cálculos de energia dos métodos estatísticos mais básicos ( t test , ANOVA etc.) podem ser concluídos em algum software (como jamovi ). No entanto, como método estatístico avançado, o cálculo de potência do modelo misto linear (LMM) não pode ser realizado em software como jamovi com interface visual e pacotes de dados e funções relacionados em R ainda são necessários. Este artigo apresenta principalmente como executar o cálculo de potência do LMM em R , incluindo: Introdução aos dados e suas funções, cenários aplicáveis e restritos e demonstrações de exemplo específicas.
simr é um pacote mais conveniente para a análise de energia LMM em R Na análise, os dados do modelo são simulados pela simulação de Monte Carlo e a proporção significativa é calculada no número de simulações definidas. A função principal é powerSim() .
Calcule o poder do efeito posterior, ou seja, calcule a potência com base nos efeitos no modelo construído com base nos dados existentes;
Calcule o poder do efeito a priori, ou seja, calcule o poder com base no tamanho de um efeito na pesquisa existente;
Estime o número de sujeitos/itens, porque o tamanho da energia é afetado pelo número de sujeitos ou itens. O modelo específico é que, quanto mais número, maior será a potência, o que nos fornece a função de estimar o número de sujeitos/itens necessários para experimentos formais com base nos efeitos do pré-experiência;
Aqui, usamos os dados de um experimento de psicologia (endereço de download de dados) para demonstrar. Este experimento coletou as reações de 10 indivíduos ( subj ) ao executar uma tarefa cognitiva ( DV ). O experimento teve duas condições ( CondA e CondB ), que foram projetados para 2*2 sujeitos e dentro do projeto ( item ). Ou seja, dois fatores aleatórios são considerados ao mesmo tempo.
Como nos concentramos em efeitos fixos e efeitos não aleatórios aqui, consideramos apenas interceptações aleatórias na parte de efeitos aleatórios ao modelar, e a parte de efeitos fixos examina os dois principais efeitos e suas interações.
A modelagem é a seguinte:
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item))
A parte do efeito fixo do modelo é:
> 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
Aqui, tomamos o efeito fixo de examinar a interação entre os dois ( CondAA2:CondBB2 ) como exemplo:
> 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
O poder do efeito fixo da interação foi 0.46 (intervalo de confiança 0.318 – 0.607 ).
Da mesma forma, o poder que calcula o efeito fixo do CondA pode ser implementado através do seguinte código:
> PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
A parte principal do efeito/interação do modelo é:
> 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
Aqui fazemos examinando o efeito principal de CondA como exemplo:
> 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
O poder do efeito principal do CondA é 0.12 (o intervalo de confiança é 0.045 – 0.243 ).
Se um tamanho de efeito relativamente estável tiver sido obtido em estudos anteriores, por exemplo, é constatado que o efeito fixo do CondA é geralmente em torno de 50, a energia também pode ser calculada com base nos efeitos anteriores. As etapas são as seguintes:
> 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
Neste momento, o poder de efeito fixo do 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
O poder principal do poder CondA é 0.94 .
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item)) # 这里仍以后验power为例
O poder do efeito fixo da interação é 0.46 e queremos testar o número de sujeitos, a potência pode atingir 0.8 . As etapas são as seguintes:
> 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
Neste momento, a energia aumentou para 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
Observe que o valor na frente de cada poder é impreciso no momento, e o valor real é:
> Pcurve$nlevels
# [1] 3 6 9 12 15 18 21 24 27 30
Pode -se observar que, quando cerca de 24 indivíduos atingem a energia, a energia pode atingir mais de 0.8 .
A faixa de tamanho da amostra pode ser definida manualmente:
> 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
O cálculo de potência do efeito no modelo generalizado é semelhante ao acima, as principais diferenças são:
method em fixed() deve ser definido como 'z' porque o efeito fixo no modelo generalizado é o teste z;method em fixed() deve ser definido como 'chisq' , porque o principal efeito no modelo generalizado é um teste qui-quadrado
