Power analysis of HLM LMM
1.0.0
Como la reproducibilidad de la psicología es cada vez más valorada, la demanda de cálculos de energía para el análisis estadístico también está aumentando. Los cálculos de potencia de los métodos estadísticos más básicos ( t test , ANOVA , etc.) se pueden completar en algún software (como jamovi ). Sin embargo, como un método estadístico avanzado, el cálculo de potencia del modelo mixto lineal (LMM) no se puede realizar en un software como jamovi con interfaz visual, y todavía se requieren paquetes de datos relacionados y funciones en R Este artículo presenta principalmente cómo realizar el cálculo de potencia de LMM en R , que incluye: Introducción a los datos y sus funciones, escenarios aplicables y restringidos, y demostraciones de ejemplo específicas.
simr es un paquete que es más conveniente para el análisis de potencia LMM en R En el análisis, los datos del modelo se simulan por simulación de Monte Carlo y la proporción significativa se calcula en el número establecido de simulaciones. La función principal es powerSim() .
Calcule la potencia del efecto posterior, es decir, calcule la potencia en función de los efectos en el modelo construido en función de los datos existentes;
Calcule la potencia del efecto a priori, es decir, calcule la potencia en función del tamaño de un efecto en la investigación existente;
Estime el número de sujetos/elementos, porque el tamaño de la potencia se ve afectado por el número de sujetos o elementos. El modelo específico es que cuanto más número es el poder, lo que nos proporciona la función de estimar el número de sujetos/ítems requeridos para experimentos formales en función de los efectos del pre-experiencia;
Aquí usamos los datos de un experimento de psicología (dirección de descarga de datos) para demostrar. Este experimento recopiló las reacciones de 10 sujetos ( subj ) cuando realiza una tarea cognitiva ( DV ). El experimento tuvo dos condiciones ( CondA y CondB ), que fueron diseñadas para 2*2 sujetos y dentro del proyecto ( item ). Es decir, dos factores aleatorios se consideran al mismo tiempo.
Debido a que nos centramos en los efectos fijos y los efectos no aleatorios aquí, solo consideramos intercepciones aleatorias en la parte de los efectos aleatorios al modelar, y la parte de efectos fijos examina los dos efectos principales y sus interacciones.
El modelado es el siguiente:
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item))
La parte de efecto fijo del modelo es:
> 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
Aquí tomamos el efecto fijo de examinar la interacción entre los dos ( CondAA2:CondBB2 ) como ejemplo:
> 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 potencia del efecto fijo de la interacción fue 0.46 (intervalo de confianza 0.318 – 0.607 ).
Del mismo modo, la potencia que calcula el efecto fijo de CondA se puede implementar a través del siguiente código:
> PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
El efecto principal/parte de interacción del modelo es:
> 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
Aquí tomamos examinando el efecto principal de CondA como ejemplo:
> 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
El poder del efecto principal de CondA es 0.12 (el intervalo de confianza es 0.045 – 0.243 ).
Si se ha obtenido un tamaño de efecto relativamente estable en estudios anteriores, por ejemplo, se encuentra que el efecto fijo de CondA generalmente es de alrededor de 50, la potencia también se puede calcular en función de los efectos anteriores. Los pasos son los siguientes:
> 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
En este momento, la potencia de efecto fijo de CondA es 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
El poder principal de la potencia de CondA es 0.94 .
> Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item)) # 这里仍以后验power为例
El poder del efecto fijo de la interacción es 0.46 , y queremos probar el número de sujetos, la potencia puede alcanzar 0.8 . Los pasos son los siguientes:
> 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
En este momento, el poder ha aumentado a 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
Tenga en cuenta que el valor frente a cada potencia es inexacto en este momento, y el valor real es:
> Pcurve$nlevels
# [1] 3 6 9 12 15 18 21 24 27 30
Se puede ver que cuando alrededor de 24 sujetos alcanzan la potencia, la potencia puede alcanzar más de 0.8 .
El rango de tamaño de la muestra se puede configurar 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
El cálculo de potencia del efecto en el modelo generalizado es similar al anterior, las principales diferencias son:
method en fixed() debe establecerse en 'z' porque el efecto fijo en el modelo generalizado es la prueba z;method en fixed() debe establecerse en 'chisq' , porque el efecto principal en el modelo generalizado es una prueba de chi-cuadrado
