效应量的计算——Cohen's d statistic
效应量的计算——Cohen's d statistic
上文ISME-人类微生物多样性与疾病的关系中提到了，采用Cohen's d statistic对效应量进行了检验。本文对此进行解释。
在统计学中，效应量（effect size）是对现象量级的定量度量。其包括两个变量之间的相关性，回归中的回归系数，平均差，甚至是发生某事的风险：如有多少人在心脏病发作后幸存下来。
效应量也是用抽样误差来估计的。大多数情况下效应量的绝对值越大，效应就越强，但是优势率（odds ratio）这个指标除外。效应量是对统计检验的补充，可以评价统计的强度。
特别是在meta分析中，其目的是将多个效应量结合起来，因此效应量的标准误差(S.E.)至关重要。
效果量的结果有助于解释研究的实质意义，而不是统计意义的显著性。
与各种统计检验的关系
基于样本的效应量不同于假设检验中使用的检验统计量，因为它们只是估计了关系的强度，而没有指定一个显著性水平来反映所观察到的关系的大小是否可能是偶然的。效果量并不直接决定显著性水平。 如果样本量足够大，非空统计比较总是会得到统计上显著的结果，除非总体效应量恰好为零(即使效应量恰好为零，也会由第一类错误而显示统计上的显著性)。例如样本量为1000，皮尔逊相关系数为0.01的样本也会具有统计学意义。
标准的效应量包括r, Cohen's d和odds ratio。不标准的效应量包括组间均值或非标准化回归系数的差异。以下情况要使用标准的效应量： 1.研究变量的度量没有内在意义 2.综合多项研究 3.综合不同尺度的研究 4.群落中由于变化而产生的效应大小
目前已经有50~100种表征效应量的方法； 如皮尔森相关的r值，r = 0.1，0.3，0.5分别被定义为small,medium,large的效应量。 本文主要关注的Cohen's d，d = 0.01,0.2,0.5,0.8,1.2,2.0分别被定义为very small,small,medium,large,very large和huge的效应量。
R中计算Cohen's d-statistic
#安装
>install.packages("compute.es")
>library(compute.es)
#主要函数tes(),对t检验的结果进行转化，得到效应量的平均差异d,d的无偏估计g, 相关系数r，Fisher's z，和对数差异比（log odds ratio）。
#还计算了这些值的方差、置信区间和p值。
#另外还有NNT (number needed to treat), U3 (Cohen's U_(3) overlapping proportions of distributions), CLES (Common Language Effect Size) and Cliff's Delta.
#已知一个t检验的t值为1.74，处理组样本30个，对照组样本31个：
>tes(t=1.74, n.1=30, n.2=31)
Mean Differences ES:
d [ 95 %CI] = 0.45 [ -0.07 , 0.96 ]
var(d) = 0.07
p-value(d) = 0.09
U3(d) = 67.21 %
CLES(d) = 62.37 %
Cliff's Delta = 0.25
g [ 95 %CI] = 0.44 [ -0.07 , 0.95 ]
var(g) = 0.07
p-value(g) = 0.09
U3(g) = 67 %
CLES(g) = 62.21 %
Correlation ES:
r [ 95 %CI] = 0.22 [ -0.04 , 0.45 ]
var(r) = 0.02
p-value(r) = 0.09
z [ 95 %CI] = 0.22 [ -0.04 , 0.49 ]
var(z) = 0.02
p-value(z) = 0.09
Odds Ratio ES:
OR [ 95 %CI] = 2.24 [ 0.88 , 5.75 ]
p-value(OR) = 0.09
Log OR [ 95 %CI] = 0.81 [ -0.13 , 1.75 ]
var(lOR) = 0.22
p-value(Log OR) = 0.09
Other:
NNT = 6.85
Total N = 61
#更加复杂一点，假设存在5个t检验：
>dat <- data.frame(id=1:5,t=rnorm(5, 2, .5),
n.t=round(rnorm(5, 25),0),
n.c=round(rnorm(5, 25),0))
id t n.t n.c
1 1 2.318254 25 24
2 2 1.943387 27 26
3 3 1.708333 25 26
4 4 1.631865 25 26
5 5 1.795588 24 24
#tes分析，结果会以数据框形式返回
>tes(t=t, n.1=n.t, n.2=n.c, level=95, dig=2, id=id, data=dat)
id N.total n.1 n.2 d var.d l.d u.d U3.d cl.d cliffs.d pval.d g var.g
1 1 49 25 24 0.66 0.09 0.07 1.25 74.62 68.03 0.36 0.03 0.65 0.08
2 2 53 27 26 0.53 0.08 -0.03 1.10 70.33 64.71 0.29 0.06 0.53 0.08
3 3 51 25 26 0.48 0.08 -0.09 1.05 68.39 63.25 0.26 0.10 0.47 0.08
4 4 51 25 26 0.46 0.08 -0.11 1.03 67.62 62.67 0.25 0.11 0.45 0.08
5 5 48 24 24 0.52 0.09 -0.07 1.11 69.79 64.30 0.29 0.08 0.51 0.08
l.g u.g U3.g cl.g pval.g r var.r l.r u.r pval.r fisher.z var.z l.z
1 0.07 1.23 74.28 67.76 0.03 0.32 0.02 0.04 0.56 0.03 0.33 0.02 0.04
2 -0.03 1.08 70.06 64.51 0.06 0.26 0.02 -0.02 0.50 0.06 0.27 0.02 -0.02
3 -0.09 1.03 68.12 63.05 0.10 0.24 0.02 -0.05 0.49 0.10 0.24 0.02 -0.05
4 -0.11 1.01 67.37 62.49 0.11 0.23 0.02 -0.06 0.48 0.12 0.23 0.02 -0.06
5 -0.07 1.09 69.49 64.08 0.08 0.26 0.02 -0.04 0.51 0.09 0.26 0.02 -0.04
u.z OR l.or u.or pval.or lOR l.lor u.lor pval.lor NNT
1 0.63 3.33 1.14 9.70 0.03 1.20 0.13 2.27 0.03 4.37
2 0.55 2.63 0.95 7.29 0.06 0.97 -0.05 1.99 0.06 5.58
3 0.53 2.38 0.85 6.71 0.10 0.87 -0.17 1.90 0.10 6.32
4 0.52 2.29 0.81 6.45 0.11 0.83 -0.21 1.86 0.11 6.65
5 0.56 2.56 0.88 7.48 0.08 0.94 -0.13 2.01 0.08 5.77还发现了另外几个包也可以做相同的分析：
####直接用函数计算
set.seed(45) ## be reproducible
>x <- rnorm(10, mean=10, sd = 1)
>y <- rnorm(10, mean = 5, sd = 5)
>cohens_d <- function(x, y) {
lx <- length(x)- 1
ly <- length(y)- 1
md <- abs(mean(x) - mean(y)) ## mean difference (numerator)
csd <- lx * var(x) + ly * var(y)
csd <- csd/(lx + ly)
csd <- sqrt(csd) ## common sd computation
cd <- md/csd ## cohen's d
}
>res <- cohens_d(x, y)
>res
[1] 0.5199662
###lsr这个包也可以
>install.packages("lsr")
>library(lsr);?lsr
>set.seed(45)
>x <- rnorm(10, 10, 1)
>y <- rnorm(10, 5, 5)
>cohensD(x,y)
[1] 0.5199662
###sffsize也可以
>install.packages("effsize")
>library(effsize); ?effsize
>set.seed(45)
>x <- rnorm(10, 10, 1)
>y <- rnorm(10, 5, 5)
>cohen.d(x,y)
Cohen's d
d estimate: 0.5199662 (medium)
95 percent confidence interval:
inf sup
-0.4353393 1.4752717
###pwr也可以,略过Reference https://en.wikipedia.org/wiki/Effect_size#Cohen's_d https://stackoverflow.com/questions/15436702/estimate-cohens-d-for-effect-size