上次介绍了两组发生率
的样本量
计算方法,通过pwr
包进行计算非常简单,可以有效地减少我们的工作量。😘
有时候我们想比较两组之间的均值
,如何计算样本量
又一次成了老大难问题。🤒
本期我们还是基于pwr
包,试一下通过两组的均值
进行样本量的估算。😏
rm(list = ls())
library(pwr)
library(tidyverse)
还是假设我们正在进行一项RCT
研究,旨在评估Treatment A
和Treatment B
之间血红蛋白A1c (HbA1c
)相对于基线
的平均变化
的差异。🤪
我们先提出研究假设,H_0 和H_1 :👇
Treatment A
和Treatment B
间HbA1c
相对于基线的平均变化没有差异。Treatment A
和Treatment B
间HbA1c
相对于基线的平均变化存在差异。接着我们还有几个参数需要设置:👇
alpha level
(通常为two-sided
);effect size
(Cohen’s d
);power
(通常为80%
);Treatment A
和Treatment B
的HbA1c
平均变化我们还是需要检索既往文献,然后做出假设。😉
这里我们假设Treatment A
的预期平均变化为1.5%
,标准差为0.25%
,Treatment B
的预期平均变化为1.0%
,标准差为0.20
。🥳
首先我们计算一下并合标准偏差
(pooled standard deviation
, σpooled
)。😂
sd1 <- 0.25
sd2 <- 0.30
sd_pooled <- sqrt((sd1^2 +sd2^2) / 2)
sd_pooled
得到了σpooled
,我们就可以计算Cohen’s d
了。🥰
mu1 <- 1.5
mu2 <- 1.0
d <- (mu1 - mu2) / sd_pooled
d
现在,我们可以利用pwr
包计Treatment A
和Treatment B
之间平均HbA1c
变化差0.5%
(1.5%
- 1.0%
)所需的样本量
,具有 80%
的power
和0.05
的显著性
。
n_i <- pwr.t.test(d = d, power = 0.80, sig.level = 0.05)
n_i
Note! 这里我们得到每组需要6个
受试者。🤔
接着是效力分析
(Power Analysis
),还是和之前的一样,如果结果不好的话,我们应该修改
或者直接终止
实验。🫠
我们绘制一下Power
随着样本量
的变化吧。🥳
n <- seq(1, 10, 1)
nchange <- pwr.t.test(d = d, n = n, sig.level = 0.05)
nchange.df <- data.frame(n, power = nchange$power * 100)
nchange.df
plot(nchange.df$n,
nchange.df$power,
type = "b",
xlab = "Sample size, n",
ylab = "Power (%)")
随着样本量的增加,power
也是不断增加的。😗
与我们之前的例子一样,随着我们增加样本量,我们减少了估算的不确定性
。🤤
通过减少这种不确定性
,有效减少了犯II类错误的可能性。🥸
接着我们固定一下两组各位6例
受试者,alpha
为0.05
。😬
我们再改变一下Treatment A
的HbA1c
相对于基线的平均变化设为0-2%
,间隔为0.1%
。🤒
mu1 <- seq(0.0, 2.0, 0.1)
d <- (mu1 - mu2) / sd_pooled
power1 <- pwr.t.test(d = d, n = 6, sig.level = 0.05)
powerchange <- data.frame(d, power = power1$power * 100)
powerchange
随着Cohen's d
的绝对值增大,power
也随之增大。⤴️
也就是说在固定样本大小的情况小,两组间差异越大我们越有可能得到阳性结果。🤠
相反,如果差异很小,那么我们在当前样本量6
的情况下没有足够的把握得到阳性结果。🤒
plot(powerchange$d,
powerchange$power,
type = "b",
xlab = "Cohen's d",
ylab = "Power (%)")
最后祝大家早日不卷!~