r - 使用 powerCurve(simr 包)进行功耗分析会产生令人困惑的输出

Question

在以下示例中，我对以下数据集执行功效分析:

hh <- data.frame(Species=c(rep("SpA", 7),rep("SpB", 5),rep("SpC", 14),rep("SpD", 10),rep("SpE", 1)),
    Skull.length=c(13.100, 14.700, 14.200, 15.400, 15.300, 15.100, 15.200, 11.100, 11.500, 12.900, 12.500, 12.400, 12.700, 12.100, 13.200, 12.300, 11.335, 12.900, 12.500, 13.190, 12.900, 14.400, 14.400, 14.300, 14.100, 14.300, 12.600, 12.900, 12.900, 14.260, 13.670, 14.720, 14.440, 14.440, 15.350, 14.970, 10.300),
    Spine.length=c(59.200, 60.100, 60.600, 67.010, 70.000, 70.300, 70.800, 53.300, 53.800, 54.200, 54.300, 56.900, 55.300, 56.600, 57.800, 57.800, 58.365, 59.900, 60.000, 60.100, 60.200, 62.900, 63.600, 63.700, 66.200, 66.700, 55.300, 55.500, 59.300, 59.740, 61.330, 65.400, 65.600, 65.800, 66.650, 68.030, 52.100))

我需要这些包:

library(lme4)
library(lmerTest) # a pimped-up version of lme4 which also provides pseudo-p-values.
library(MuMIn) # gives pseudo-R-squared via r.squaredGLMM()
library(pwr) # power analysis for lm
library(simr) # power analysis for generalized linear mixed models by simulation 

如果我要测试 Skull.length 之间的相关性和Spine.length忽略Species的作用我会这样做:

lm1 <- lm(Skull.length~Spine.length, data=hh)
summary(lm1)$adj.r.squared # 0.7696584

然后使用包 pwr 进行功率分析来测试我的样本量是否足够大将很容易:

p.out <- pwr.r.test(r = sqrt(summary(lm1)$adj.r.squared), sig.level = 0.05, power = 0.8, alternative = "greater")
# To detect r = 0.8773018 or greater with sig.level = 0.05 and power = 0.8, n >= 6 is required

但我想考虑hh$Species如下图所示:

mem.skull.vs.body <- glmer(Skull.length ~ Spine.length + (1| Species),
                            data=hh,
                            family="gaussian")

产生:

Fixed effects:
             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)   0.73958    1.32239 23.50147   0.559    0.581    
Spine.length  0.20848    0.02173 22.72726   9.593 1.87e-09 ***

[模型参数的数据和线性回归 mem.skull.vs.body ]

我的模型的斜率，0.20848 ，是我对效应大小的衡量。要找出检测至少 0.1 效应量所需的样本量:

fixef(mem.skull.vs.body)["Spine.length"] <- 0.1
powerSim(mem.skull.vs.body, nsim=1000)

这给出:

Power for predictor 'Spine.length', (95% confidence interval):
      98.90% (98.04, 99.45)

这表明我的样本量(37 个人，每个人来自五个物种之一)对于我正在测试的模型来说足够了，但是当我继续使用 powerCurve(mem.skull.vs.body, nsim=1000) 进行仔细检查时我得到:

Power for predictor 'Spine.length', (95% confidence interval),
by largest value of Spine.length:
   53.8:  0.00% ( 0.00,  0.37) - 3 rows
   55.3:  5.40% ( 4.08,  6.99) - 7 rows
   57.8:  5.20% ( 3.91,  6.76) - 12 rows
   59.3: 12.30% (10.33, 14.50) - 15 rows
   60.1: 21.50% (18.99, 24.18) - 20 rows
  61.33: 30.60% (27.75, 33.56) - 23 rows
   65.4: 61.40% (58.30, 64.43) - 27 rows
   66.2: 80.00% (77.38, 82.44) - 30 rows
  68.03: 94.80% (93.24, 96.09) - 34 rows
   70.8: 98.40% (97.41, 99.08) - 37 rows

下面是上述值的图表:

我发现这个输出即使不可疑也令人困惑，因为:

• 这表明我需要超过 65 个观察值的样本才能获得 与powerSim()中的估计相比，检测到效应大小为 0.1 的可能性为 80% ;
• x 轴的值范围非常接近 hh$Spine.length 假设的值范围，介于 52.1 和 70.8 之间。

它看起来非常像函数 powerCurve在其默认设置中，会将 x 值的大小与样本大小混淆。有没有办法更改 powerCurve 的设置以避免这种困惑？

---

更新(2019 年 4 月):

自从我提出这个问题以来，软件包开发人员已经修改了函数 powerCurve以反射(reflect) pete 下面提供的解释.

Answer 1

powerCurve 采用 along 参数，默认为第一个固定协变量。并非所有变量都有意义，如本例所示。

在这种情况下，您可以添加一个“观察”变量并沿着该变量运行功效曲线:

hh$obs <- 1:37
pc <- powerCurve(mem.skull.vs.body, along="obs")

然后plot(pc)会给出更直观的结果。

---

如果您想更多地控制绘图，我建议使用summary来获取原始数字，然后根据您认为合适的方式绘制它们。请注意，nrow 列当前仅在 github 版本中可用(如果您将来阅读本文，则在版本 > 1.0.5 中可用)。

summary(pc)
#    nrow nlevels successes trials mean     lower      upper
# 1     3       3         0    100 0.00 0.0000000 0.03621669
# 2     7       7         0    100 0.00 0.0000000 0.03621669
# 3    11      11         9    100 0.09 0.0419836 0.16398226
# 4    14      14        18    100 0.18 0.1103112 0.26947709
# 5    18      18        32    100 0.32 0.2302199 0.42076686
# 6    22      22        67    100 0.67 0.5688272 0.76080147
# 7    26      26        90    100 0.90 0.8237774 0.95099531
# 8    29      29        91    100 0.91 0.8360177 0.95801640
# 9    33      33        98    100 0.98 0.9296161 0.99756866
# 10   37      37        98    100 0.98 0.9296161 0.99756866