模型可视化>> Power-analysis-of-HLM-LMM>> 返回
项目作者： usplos

项目描述 ：
Estimate the power of linear mixed model
高级语言：
项目地址: git://github.com/usplos/Power-analysis-of-HLM-LMM.git
创建时间： 2019-11-09T23:50:34Z
项目社区：https://github.com/usplos/Power-analysis-of-HLM-LMM
开源协议：
下载


线性混合模型的power分析

随着心理学的可重复性日渐受到重视，对统计分析的power计算的需求也日益增多，
多数基础统计方法 (t test, ANOVA 等) 的power计算可以在某些软件 (比如jamovi) 中完成。
但作为一种高级统计方法，线性混合模型 (linear mixed model, LMM) 的power计算尚不能在jamovi这样具有可视化界面的软件中进行，
仍需利用R中的相关数据包及函数。本文主要介绍在 R 中如何进行LMM的power计算，内容包括：数据及其函数介绍，可适用及受限制的情景，
以及具体的示例演示。

数据包及其函数

simr 数据包是 R 中进行LMM的power分析较为便利的包。
在分析中，依靠蒙泰卡罗模拟的方式来模拟模型数据，并且计算在设定的模拟次数中，达到显著的比例。主要的函数为powerSim()。

适用情景

  • 计算后验效应的power，即根据基于已有数据所建模型中的效应来计算power；

  • 计算先验效应的power，即根据已有研究中某效应的大小来计算power；

  • 预估被试量/项目量，因为power的大小受到被试或项目数量的影响，具体模式为数量越多，power越大，
    这为我们提供了根据预实验的效应来估计正式实验所需被试量/项目量的功能；

限制情景

  • 整个power分析耗时较长（目前该数据包作者也在试图解决这个问题，用户也可以采用多线程计算的方式来缓解）
  • 即使先验效应的power，也需要预实验的数据

示例

计算后验效应的power

这里我们利用一个心理学实验的数据(数据下载地址)来演示，
该实验收集了10名被试(subj)在执行一项认知任务时的反应时(DV)，实验有两个条件(CondA和CondB)，为2*2被试内和项目(item)内设计。
即同时考虑被试和项目两个随机因子。

因为这里我们关注固定效应，非随机效应，因此建模时在随机效应部分只考虑随机截距，固定效应部分考察两个主效应及其交互作用，

建模如下：

  1. > Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item))

考察固定效应的power

模型的固定效应部分为:

  1. > summary(Model)$coef
  2. #                  Estimate Std. Error        df    t value     Pr(>|t|)
  3. # (Intercept)     279.43090   23.37537  11.71977 11.9540740 6.422210e-08
  4. # CondAA2          24.29565   13.49694 498.43142  1.8000866 7.245150e-02
  5. # CondBB2          12.18484   13.36445 509.31328  0.9117351 3.623395e-01
  6. # CondAA2:CondBB2 -32.82881   19.29629 509.76684 -1.7013020 8.949610e-02

这里以考察二者交互作用的固定效应(CondAA2:CondBB2)为例：

  1. > PowerAB_ttest = simr::powerSim(fit = Model, # 要考察的模型
  2.                                  test = fixed('CondAA2:CondBB2', # 要考察的固定效应的名称
  3.                                               method = 't'), # 选取检验方法，因为固定效应为t检验，因此method设置为t
  4.                                  nsim=50) # 设置模拟次数，建议设置为500 (此时可以获取到较稳定的power)
  1. > PowerAB_ttest
  2. # Power for predictor 'CondAA2:CondBB2', (95% confidence interval):
  3. #       46.00% (31.81, 60.68)
  4. # 
  5. # Test: t-test with Satterthwaite degrees of freedom (package lmerTest)
  6. #       Effect size for CondAA2:CondBB2 is -33.
  7. # 
  8. # Based on 50 simulations, (2 warnings, 0 errors)
  9. # alpha = 0.05, nrow = 553
  10. # 
  11. # Time elapsed: 0 h 0 m 8 s
  12. # 
  13. # nb: result might be an observed power calculation

交互作用固定效应的power为0.46(置信区间为0.318 – 0.607)。

同理，计算CondA的固定效应的power可通过以下代码实现：

  1. > PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)

计算主效应/交互作用的power

模型的主效应/交互作用部分为：

  1. > anova(Model)
  2. # Type III Analysis of Variance Table with Satterthwaite's method
  3. #             Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
  4. # CondA         8546    8546     1 489.18  0.6694 0.4137  
  5. # CondB         2453    2453     1 509.78  0.1921 0.6613  
  6. # CondA:CondB  36951   36951     1 509.77  2.8944 0.0895 .
  7. # ---
  8. # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这里以考察CondA主效应为例：

  1. > PowerA_Ftest = simr::powerSim(fit = Model, # 要考察的模型
  2.                                 test = fixed('CondA', # 要考察的主效应/交互作用的名称
  3.                                              method = 'f'), # 选取检验方法，因为主效应为f检验，因此method设置为f
  4.                                 nsim=50) # 设置模拟次数，建议设置为500 (此时可以获取到较稳定的power)
  1. > PowerA_Ftest
  2. # Power for predictor 'CondA', (95% confidence interval):
  3. #       12.00% ( 4.53, 24.31)
  4. # 
  5. # Test: Type-II F-test (package car)
  6. # 
  7. # Based on 50 simulations, (50 warnings, 0 errors)
  8. # alpha = 0.05, nrow = 553
  9. # 
  10. # Time elapsed: 0 h 0 m 58 s
  11. # 
  12. # nb: result might be an observed power calculation

CondA的主效应的power为0.12(置信区间为0.045 – 0.243)。

计算先验效应的power

如果先前研究中已经得到了较稳定的效应量，比如发现CondA的固定效应一般在50左右,也可以根据先验的效应来计算power。步骤如下：

  1. 将模型的后验固定效应改为先验固定效应
  1. > fixef(Model) # 查看模型固定效应
  2. #     (Intercept)         CondAA2         CondBB2 CondAA2:CondBB2 
  3. #       279.43090        24.29565        12.18484       -32.82881
  1. > fixef(Model)[2] = 50 # 修改CondA的效应
  2. > fixef(Model) # 查看修改后的固定效应
  3. #     (Intercept)         CondAA2         CondBB2 CondAA2:CondBB2 
  4. #       279.43090        50.00000        12.18484       -32.82881
  1. 根据先验固定效应计算power(与计算后验power的方法相同)
  1. > PowerA_ttest = simr::powerSim(fit = Model, test = fixed('CondAA2', method = 't'), nsim=50)
  2. > PowerA_ttest
  3. # Power for predictor 'CondAA2', (95% confidence interval):
  4. #       88.00% (75.69, 95.47)
  5. # 
  6. # Test: t-test with Satterthwaite degrees of freedom (package lmerTest)
  7. #       Effect size for CondAA2 is 50.
  8. # 
  9. # Based on 50 simulations, (3 warnings, 0 errors)
  10. # alpha = 0.05, nrow = 553
  11. # 
  12. # Time elapsed: 0 h 0 m 9 s

此时CondA的固定效应power为0.88.

  1. > PowerA_Ftest = simr::powerSim(fit = Model, test = fixed('CondA', method = 'f'), nsim=50)
  2. > PowerA_Ftest
  3. # Power for predictor 'CondA', (95% confidence interval):
  4. #       94.00% (83.45, 98.75)
  5. # 
  6. # Test: Type-II F-test (package car)
  7. # 
  8. # Based on 50 simulations, (50 warnings, 0 errors)
  9. # alpha = 0.05, nrow = 553
  10. # 
  11. # Time elapsed: 0 h 1 m 4 s

CondA的主效应power为0.94.

预估被试量

  1. > Model = lmer(data = DemoData, DV~CondA*CondB+(1|subj)+(1|item)) # 这里仍以后验power为例

前面已知交互作用固定效应的power为0.46，我们想考察被试量为多少时，该power可以达到0.8。步骤如下：

  1. 建立扩充被试后的模型，比如这里扩充为30名被试
  1. > Model2 = extend(object = Model, along = 'subj', n = 30)
  1. 计算此时的power
    ```

    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

```

此时power已经升到0.84。

  1. 如果想知道具体的power随被试量增大的变化趋势，可利用powerCurve()函数
  1. > Pcurve = powerCurve(fit = Model2, test = fixed('CondAA2', method = 't'), along = 'subj', nsim=50)
  3. > Pcurve
  4. # Power for predictor 'CondAA2', (95% confidence interval), by largest value of subj:
  5. #      11: 22.00% (11.53, 35.96) - 180 rows
  6. #      14: 36.00% (22.92, 50.81) - 325 rows
  7. #      17: 50.00% (35.53, 64.47) - 507 rows
  8. #       2: 58.00% (43.21, 71.81) - 667 rows
  9. #      22: 64.00% (49.19, 77.08) - 837 rows
  10. #      25: 58.00% (43.21, 71.81) - 989 rows
  11. #      28: 74.00% (59.66, 85.37) - 1166 rows
  12. #      30: 82.00% (68.56, 91.42) - 1318 rows
  13. #       6: 86.00% (73.26, 94.18) - 1489 rows
  14. #       9: 90.00% (78.19, 96.67) - 1659 rows
  15. # 
  16. # Time elapsed: 0 h 1 m 39 s

注意，此时每个power前面的数值不准确，真实的数值为:

  1. > Pcurve$nlevels
  2. # [1]  3  6  9 12 15 18 21 24 27 30

可以看到大概到24名被试时，power可以到0.8以上.

可以人为设置样本量的范围:

  1. > Pcurve = powerCurve(fit = Model2, test = fixed('CondAA2', method = 't'), along = 'subj', nsim=50, 
  2.                       breaks=c(21, 22, 23, 24, 25, 30)) # 考察样本量为21、22、23、24、25、30时的power
  4. > Pcurve
  5. # Power for predictor 'CondAA2', (95% confidence interval), by largest value of subj:
  6. #      28: 72.00% (57.51, 83.77) - 1166 rows
  7. #      29: 74.00% (59.66, 85.37) - 1220 rows
  8. #       3: 70.00% (55.39, 82.14) - 1262 rows
  9. #      30: 74.00% (59.66, 85.37) - 1318 rows
  10. #       4: 76.00% (61.83, 86.94) - 1369 rows
  11. #       9: 80.00% (66.28, 89.97) - 1659 rows
  12. # 
  13. # Time elapsed: 0 h 1 m 12 s
  15. > Pcurve$nlevels
  16. # [1] 21 22 23 24 25 30

计算广义线性混合模型中效应的power

广义模型中效应的power计算与上面的类似，主要区别是：

  • 考察固定效应时，fixed()中method参数应设置为'z'，因为广义模型中固定效应为z检验;
  • 考察主效应时，fixed()中method参数应设置为'chisq'，因为广义模型中主效应为卡方检验