1写在前面
估算前瞻性研究
的样本量
是我们在招募受试者之前首先要做的事情之一。😘
招募受试者太少会无法得到准确的答案,招募的太多又是巨大的浪费,所以需要估算最佳的受试者数量。🧐
本期我们介绍一下如何使用pwr
包进行样本量的估算。😏
2用到的包
rm(list = ls())
library(pwr)
library(tidyverse)
3研究假设
假设我们准备进行一个RCT
研究,研究Treatment A
和Treatment B
的疗效,结局事件为Response
或No response
的二分类
结局。🤫
那我们现在就有了研究假设,H_0 和H_1 了:👇
- H_0:
Treatment A
和Treatment B
间结局事件无差异。 - H_1 :
Treatment A
和Treatment B
间结局事件有差异。
通常我们还有几个参数需要设置:👇
alpha level
(通常为two-sided
);effect size
(h
);power
(通常为80%
)
4计算样本量
这里我们需要用到pwr.2p.test
函数,我们通常需要设置Treatment A
和Treatment B
的response
比例,这个大家可以通过既往的文献来查找。😘
如果你做的研究非常新,在过去的文献中找不到的话,可以假设为50%
。😂
这里我们假设Treatment A
反应率是60%
, Treatment B
反应率是50%
,这样Treatment A
和Treatment B
间的response
比例就相差了10%
,哈哈哈哈。😉
power1 <-pwr.2p.test(h = ES.h(p1 = 0.60, p2 = 0.50), sig.level = 0.05, power = .80)
power1
Note! 需要注意的是这里n
只是一个组的数量,实际需要的总数量是需要double
的。🤞
5Power Analysis
接着是效力分析
(Power Analysis
),主要是用来确定在指定显著性条件下所需要的样本量
并评估该实验设计的统计效力
。😗
通过Power Analysis
,我们也能给出在现有的样本量下该实验结论的可靠性
。🤩
如果结论的可靠性非常低,那么几乎可以认为实验是无效的
,我们应该修改
或者直接终止
实验。🫠
plot(power1)
这里我们可以看到power
随样本量增加的变化。🤓
6改变一下
这个时候我们改一下Treatment A
反应率,并且将每组的样本量固定为388
,α值
为0.05
。
p1 <- seq(0.5, 1.0, 0.05)
power1 <-pwr.2p.test(h = ES.h(p1 = p1, p2 = 0.50),
n = 388,
sig.level = 0.05)
power1
可视化一下可以看到Treatment A
反应比例越高,和Treatment B
差异越大,power
就越大。😏
powerchange <- data.frame(p1, power = power1$power * 100)
plot(powerchange$p1,
powerchange$power,
type = "b",
xlab = "Proportion of Responders in Treatment A",
ylab = "Power (%)")
最后祝大家早日不卷!~