使用 powerCurve（package simr）的功率分析给出了令人困惑的输出
Power analysis with powerCurve (package simr) gives confusing output
在以下示例中，我对以下数据集执行功效分析：
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(); 的估计值相比，有 80% 的机会检测到 0.1 的效应大小
- x轴的数值范围非常接近
hh$Spine.length假设的数值范围，介于52.1和70.8之间。
它看起来很像函数 powerCurve 在其默认设置中混淆了 x 值的大小与样本大小。有没有办法更改 powerCurve 的设置以避免这种混淆？
更新（2019 年 4 月）：
自从我问了这个问题后，软件包开发人员修改了函数 powerCurve 以反映 pete.
下面提供的解释
powerCurve 采用默认为第一个固定协变量的 along 参数。并非所有变量都有意义，如本例所示。
在这种情况下，您可以添加一个 "observation" 变量和 运行 沿着该变量的幂曲线：
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
我想最好对 pete 的回答补充一点关于混乱的解释。
在 Marco Plebani 的模拟中，扩展是沿着“hh$Spine.length”，即“66.2”不能理解为样本量，而是脊柱的长度。
在皮特的模拟中，hh$obs 的值对应于样本数。
为了得到对应样本量的 80% 功效，我们可以稍微改进一下 pete 的解决方案：
mem.skull.vs.body2 <- update(mem.skull.vs.body, control=lmerControl(check.conv.singular = .makeCC(action = "ignore", tol = 1e-4))) #disable singular warning message
powerCurve(mem.skull.vs.body2, along="obs", breaks=c(22, 23, 24, 25,26,27))
Calculating power at 10 sample sizes along Spine.length
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: 9.90% ( 8.12, 11.92) - 7 rows
57.8: 18.90% (16.52, 21.47) - 12 rows
59.3: 48.60% (45.46, 51.75) - 15 rows
60.1: 78.30% (75.61, 80.82) - 20 rows
61.33: 92.90% (91.13, 94.41) - 23 rows
65.4: 99.50% (98.84, 99.84) - 27 rows
66.2: 100.0% (99.63, 100.0) - 30 rows
68.03: 100.0% (99.63, 100.0) - 34 rows
70.8: 100.0% (99.63, 100.0) - 37 rows
我不知道为什么我的模拟结果与皮特的相差很大。
我切换到
powerCurve(mem.skull.vs.body, along="obs", breaks=c(14,16,17,18,20,22))
Calculating power at 6 sample sizes along obs
Power for predictor 'Spine.length', (95% confidence interval),==========================================================|
by largest value of obs:
14: 47.00% (36.94, 57.24) - 14 rows
16: 61.00% (50.73, 70.60) - 16 rows
17: 83.00% (74.18, 89.77) - 17 rows
18: 98.00% (92.96, 99.76) - 18 rows
20: 100.0% (96.38, 100.0) - 20 rows
22: 100.0% (96.38, 100.0) - 22 rows
那么看来 17 个样本足以提供 >=80% 的功率。
用17个样本来验证。
library(dplyr)
hh17 <- sample_n(hh, size=17, replace=F)
model17 <- lmer(Skull.length ~ Spine.length + (1| Species), data=hh17)
powerSim(model17,nsim=100)
Power for predictor 'Spine.length', (95% confidence interval):==========================================================|
97.00% (91.48, 99.38)
Test: unknown test
Effect size for Spine.length is 0.17
Based on 100 simulations, (0 warnings, 0 errors)
alpha = 0.05, nrow = 17
Time elapsed: 0 h 0 m 9 s
上面的采样结果似乎过度优化了。通过更多的模拟，17 个样本可能足以提供超过 80% 的功效。
在以下示例中，我对以下数据集执行功效分析：
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(); 的估计值相比，有 80% 的机会检测到 0.1 的效应大小
- x轴的数值范围非常接近
hh$Spine.length假设的数值范围，介于52.1和70.8之间。
它看起来很像函数 powerCurve 在其默认设置中混淆了 x 值的大小与样本大小。有没有办法更改 powerCurve 的设置以避免这种混淆？
更新（2019 年 4 月）：
自从我问了这个问题后，软件包开发人员修改了函数 powerCurve 以反映 pete.
powerCurve 采用默认为第一个固定协变量的 along 参数。并非所有变量都有意义，如本例所示。
在这种情况下，您可以添加一个 "observation" 变量和 运行 沿着该变量的幂曲线：
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
我想最好对 pete 的回答补充一点关于混乱的解释。 在 Marco Plebani 的模拟中，扩展是沿着“hh$Spine.length”，即“66.2”不能理解为样本量，而是脊柱的长度。 在皮特的模拟中，hh$obs 的值对应于样本数。 为了得到对应样本量的 80% 功效，我们可以稍微改进一下 pete 的解决方案：
mem.skull.vs.body2 <- update(mem.skull.vs.body, control=lmerControl(check.conv.singular = .makeCC(action = "ignore", tol = 1e-4))) #disable singular warning message
powerCurve(mem.skull.vs.body2, along="obs", breaks=c(22, 23, 24, 25,26,27))
Calculating power at 10 sample sizes along Spine.length
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: 9.90% ( 8.12, 11.92) - 7 rows
57.8: 18.90% (16.52, 21.47) - 12 rows
59.3: 48.60% (45.46, 51.75) - 15 rows
60.1: 78.30% (75.61, 80.82) - 20 rows
61.33: 92.90% (91.13, 94.41) - 23 rows
65.4: 99.50% (98.84, 99.84) - 27 rows
66.2: 100.0% (99.63, 100.0) - 30 rows
68.03: 100.0% (99.63, 100.0) - 34 rows
70.8: 100.0% (99.63, 100.0) - 37 rows
我不知道为什么我的模拟结果与皮特的相差很大。 我切换到
powerCurve(mem.skull.vs.body, along="obs", breaks=c(14,16,17,18,20,22))
Calculating power at 6 sample sizes along obs
Power for predictor 'Spine.length', (95% confidence interval),==========================================================|
by largest value of obs:
14: 47.00% (36.94, 57.24) - 14 rows
16: 61.00% (50.73, 70.60) - 16 rows
17: 83.00% (74.18, 89.77) - 17 rows
18: 98.00% (92.96, 99.76) - 18 rows
20: 100.0% (96.38, 100.0) - 20 rows
22: 100.0% (96.38, 100.0) - 22 rows
那么看来 17 个样本足以提供 >=80% 的功率。 用17个样本来验证。
library(dplyr)
hh17 <- sample_n(hh, size=17, replace=F)
model17 <- lmer(Skull.length ~ Spine.length + (1| Species), data=hh17)
powerSim(model17,nsim=100)
Power for predictor 'Spine.length', (95% confidence interval):==========================================================|
97.00% (91.48, 99.38)
Test: unknown test
Effect size for Spine.length is 0.17
Based on 100 simulations, (0 warnings, 0 errors)
alpha = 0.05, nrow = 17
Time elapsed: 0 h 0 m 9 s
上面的采样结果似乎过度优化了。通过更多的模拟，17 个样本可能足以提供超过 80% 的功效。