DALS004-统计推断(Inference)02-CLT与T-test

title: DALS003-统计推断(Inference)02-CLT与T-test
date: 2019-07-22 12:0:00
type: "tags"
tags:

  • 统计推断
    categories:
  • 生物统计

中心极限定理(Central Limit Theorem)与t分布

中心极限定理(CLT)

当样本的数目非常大的，一个随机样本的均值\bar{Y}就会服从正态分布，这个正态分布的均值是总体均值(\mu_{Y})，标准差是总体的标准差\sigma_{Y}除以样本的数目N的平方根。

如果我们将一个随机变量减去一个常数，那么这个新生成的随机变量就会领衔这个常数，从数学角度来说，如果X是一个服从均值为\mu的分布，a是一个常数，那么X-a的均值的分布的均值就是\mu-a。这个也适应用乘法与标准差(SD)。

如果X是一个随机变量，它的分布的均值\mu和标准差(SD)为\sigma，而a是一个常数，那么aX这个新的随机变量分布的均值是a\mu，标准差为|a|\sigma。为了说明这个问题，我们就假设，每只小鼠的体重都减去10，那么平均体重也应该会减去10，类似的，如果我们将小鼠的体重单位由克(g, gram)换成毫克(mg, milligram)，也就相当于把原来的数字都乘以了1000，那么数字的扩散程度也会变大。

这就说明了，如果我们获得了N个样本，那么这个公式\frac{\overline{Y}-\mu}{\sigma_{Y} / \sqrt{N}}就接近于均值为0，标准差为1的正态分布。

现在我们感兴趣的是两个样本的均值的差值。这个数学公式就比较有用了。

如果我们有两个随机变量，即X和Y，它们的总体均值与总体均值差为\mu_{X}和\mu_{Y}，\sigma_{X}和\sigma_{Y}，然后我们就会得到这些结果：

X+Y的均值：\mu_{X}+\mu_{Y}。

也会得到Y-X的均值，其实这也相当于Y-X=Y+aX，其中a=-1，这其实也就是说Y-X的均值为\mu_{Y}-\mu_{X}。

如果说X和Y是相互独立的样本，那么Y+X的方差，就是\sigma_{Y}^2+\sigma_{X}^2，从这里我们可以推到，Y-X，也即Y+aX（其中a为-1）的方差是\sigma_{Y}^2+a^2\sigma_{X}^2=\sigma_{Y}^2+\sigma_{X}^2。因此说，Y-X这个差值的方差也是两个样本的方差和，这有点反直觉，我们需要知道的就是，X与Y这两个样本是相互独立的，它们确实不相互影响。

从数学角度来理解上面的内容对我们的研究目（两个样本的均值，以及均值的差值）。因此这两个样本都服从正态分布，它们的差值也服从正态分布，方差是两个样本总体的方差之和。在零假设（也就是说这两个样本所代表的两个总体的均值没有差异）成立的前提下，样本均值\bar{Y}-\bar{X}的差值的分布接近于均值为0（也就是说没有差异），标准差为\sqrt{\sigma_{X}^{2}+\sigma_{Y}^{2}} / \sqrt{N}的正态分布，也就是说下面的这个公式：
\frac{\overline{Y}-\overline{X}}{\sqrt{\frac{\sigma_{X}^{2}}{M}+\frac{\sigma_{Y}^{2}}{N}}}
的值接近于标准正态分布（均值为0，标准差为1）。

通过这种近似计算，我们就能够将计算p值的过程简单化，因为我们知道任何值在标准正态分布下的比例。例如，这些值中只有5%的比例会大于2（绝对值），计算过程如下所示：

pnorm(-2) + (1 - pnorm(2))

计算结果如下所示：

> pnorm(-2) + (1 - pnorm(2))
[1] 0.04550026

从上面结果我们就可以知道，只有5%的可能性会大于2，我们前面计算的obstdiff的值是3.020833，因此12只小鼠足够了，不需要买更多的小鼠。

但是，这还没完。因为我们并不知道总体的标准差，也就是不知道\sigma_{X}与\sigma_{Y}。它们是未知的总体参数，但是，我们可以通过样本的标准差对它们进行估算，估算的总体标准差我们不能希腊字母\mu来表示，只能用s来表示，那么它们的估计值如下所示：
s_{X}^{2}=\frac{1}{M-1} \sum_{i=1}^{M}\left(X_{i}-\overline{X}\right)^{2} \text { and } s_{Y}^{2}=\frac{1}{N-1} \sum_{i=1}^{N}\left(Y_{i}-\overline{Y}\right)^{2}
我们需要注意的是，上面的公式里面除以的是M-1与N-1，而不是前面的M与N，这有数学方面计算的原因，这里不表，记住就行。但是，为了更加直观地说明一些东西，我们还是要涉及一点。试想，如果你有2个数，把它们画在数轴上，它们之间的距离如果是L，那么中间的那一点就是它们的平均值，这个平均值离它们的距离分别都是0.5L。因此，从这个角度来说，你从这2个数中的一个就能够获取这个信息（这涉及自由度的问题），而不需要2个数，这点不怎么重要，重要的一点是，我们是用了s_{X}和s_{Y}来估计总体的均值，即\sigma_{X}和\sigma_{Y}，因为总体均值未知，只能靠样本均值来估计。

因此我们就把前面的那个公式：
\frac{\overline{Y}-\overline{X}}{\sqrt{\frac{\sigma_{X}^{2}}{M}+\frac{\sigma_{Y}^{2}}{N}}}
定义为下面的这个样子：
\frac{\overline{Y}-\overline{X}}{\sqrt{\frac{s_{X}^{2}}{M}+\frac{s_{Y}^{2}}{N}}}
如果M=N，那么公式就可以化简为如下所示：
\sqrt{N} \frac{\overline{Y}-\overline{X}}{\sqrt{s_{X}^{2}+s_{Y}^{2}}}
CLT告诉我们，当M与N足够大的时候，上面的这个随机变量就会服从均值为0，标准差为1的正态分布。

t分布

CLT的计算依赖于大量的样本，也就是我们指的渐进性结果。当CLT不适用时，还需要另外一种不依赖于渐进性结果的选择。当原始总体来源于一个变量，也就是Y时，如果这个Y是服从均值为0的标准正态分布，那么我们就可以计算下面这个变量的分布：
\sqrt{N} \frac{\overline{Y}}{s_{Y}}
这是两个随机变量的比值，因此这个比值不一定服从正态分布。事实上，这个公式的分母会偶然变小，因此就相应地会增加观察到更大值的概率。Gosset最初发现的这种分布，由于他开始在论文中以Student的名义投的稿，因此就称这种分布为Student‘s t-distribution，也就是t分布。

现在我们以小鼠的表型数据为例说明一下，我们创建2向量，一个当作对照总体，一个当作高脂总体，如下所示：

url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extd\
ata/mice_pheno.csv"
filename <- "mice_pheno.csv"
download(url,destfile=filename)
dat <- read.csv(filename)
head(dat)

library(dplyr)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% select(Bodyweight) %>% unlist

这里需要注意提，我们假设的是的y_{1},y_{2},\dots,y_{n}的分布，而不是随机变量\bar{Y}的分布。现在我们来看一下这两个总体的分布：

library(rafalib)
mypar(1,2)
hist(hfPopulation)
hist(controlPopulation)
image

我们可以使用qq图来确认一下上面的分布是比较接近正态分布的，这里我们只作初步了解，更详细的了解在后面。如果在qq图上，点大致都藻在一条直线附近，它们就属于正态分布，如下所示：

mypar(1,2)
qqnorm(hfPopulation)
qqline(hfPopulation)
qqnorm(controlPopulation)
qqline(controlPopulation)
image

练习

P47

中心极限定理的实际运用

现在我们使用上面的数据来看一下，中心极限定理是如何逼近样本的均值的：

url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extd\
ata/mice_pheno.csv"
filename <- "mice_pheno.csv"
download(url,destfile=filename)
dat <- read.csv(filename)
head(dat)

library(dplyr)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% select(Bodyweight) %>% unlist
mu_hf <- mean(hfPopulation)
mu_control <- mean(controlPopulation)
mu_hf - mu_control

结果如下所示：

> head(dat)
  Sex Diet Bodyweight
1   F   hf      31.94
2   F   hf      32.48
3   F   hf      22.82
4   F   hf      19.92
5   F   hf      32.22
6   F   hf      27.50
> 
> library(dplyr)
> controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
> hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% select(Bodyweight) %>% unlist
> mu_hf <- mean(hfPopulation)
> mu_control <- mean(controlPopulation)
> mu_hf - mu_control
[1] 2.375517

现在计算一下种群的标准差，这里不要使用R的sd()函数，因此这个函数是计算样本的标准差的，它会除以(n-1)，当我们想要对总体进行估计时，按以下代码计算：

x <- controlPopulation
N <- length(x)
populationvar <- mean((x-mean(x))^2)
populationvar
var(x)
var(x)*(N-1)/N
identical(var(x)*(N-1)/N,populationvar)
identical(var(x), populationvar)

计算结果如下所示：

> populationvar
[1] 11.67205
> var(x)
[1] 11.72416
> var(x)*(N-1)/N
[1] 11.67205
> identical(var(x)*(N-1)/N,populationvar)
[1] TRUE
> identical(var(x), populationvar)
[1] FALSE

因此，为了在数学计算上是准确的，我们不使用sd()或var()函数，我们可以使用rafalib包中的popvar()与popsd()函数，如下所示：

library(rafalib)
sd_hf <- popsd(hfPopulation)
sd_control <- popsd(controlPopulation)

这里需要注意，在实际统计中，我们并不会计算总体参数，我们通常是由样本来估计总体参数，现在我们从总体中抽样，如下所示：

N <- 12
hf <- sample(hfPopulation, N)
control <- sample(controlPopulation,N)

CLT告诉我们，对于更大的N来说，这些样本都近拉于正态分布。作为一种经验，这个N通常是大于30，这仅是一个经验法则，具体来说，它还需要取决于总体的分布。这里我们实际上可以查看一下近似值，并对各种N值进行计算，代码如下所示：

Ns <- c(3, 12, 25, 50)
B <- 10000
res <- sapply(Ns, function(n){
  replicate(B, mean(sample(hfPopulation, n)) - mean(sample(controlPopulation, n)))
})

现在使用qq图看一下CLT理论对这些数据的近似值的适应情况。如果这些数据非常接近正态分布，那么这些数据点就会落在一个直线上（这个直接是正态分布的分位数(quantiles))。越偏离，说明数据越不符合正态分布，如下所示：

mypar(2,2)
for (i in seq(along=Ns)){
  titleavg <- signif(mean(res[,i]),3)
  titlesd <- signif(popsd(res[,i]),3)
  title <- paste0("N=", Ns[i], "Avg=", titleavg, " SD=", titlesd)
  qqnorm(res[,i], main = title)
  qqline(res[,i], col=2)
  
}
image

从上面结果可以看出来，第3组数据的拟合结果最合适，这是因为它的总体相对最接近正态分布，平均值也最接近于正态分布。在实际中，我们可以计算这样一个比值：除以估计的标准误差来判断，如下所示：

Ns <- c(3, 12, 25, 50)
B <- 10000
computetstat <- function(n){
  y <- sample(hfPopulation, n)
  x <- sample(controlPopulation, n)
  (mean(y)- mean(x))/sqrt(var(y)/n+var(x)/n)
}
res <- sapply(Ns, function(n){
  replicate(B, computetstat(n))
})

mypar(2,2)
for (i in seq(along=Ns)){
  qqnorm(res[,i], main = Ns[i])
  qqline(res[,i], col=2)
  
}
image

从上面可以看出来，当N=3时，CLT并不会提供一个有用的估计，当N=12时，在较大值的点方面，稍微有点偏离。当N=25或N=50时，所有的点基本上都在直线上。

也就是说，在这个案例中，只要N=12就能证明CLT，就像前言提到的那样，在多数情况下，这种模拟并不会很好。我们在这里只是通过这种模拟手段来说明CLT背后的思想与局限。

练习

P52

t检验的实际计算

现在我们来看一下如何通过计算得到p值，先来导入数据，在这一步中，我们要确定哪些是treatment组，哪些是control组，并且计算出它们的均值差异，如下所示：

library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/femaleMiceWeights.csv")
dat <- read.csv(filename)
control <- filter(dat,Diet=="chow") %>% select(Bodyweight) %>% unlist
treatment <- filter(dat,Diet=="hf") %>% select(Bodyweight) %>% unlist
diff <- mean(treatment) - mean(control)
print(diff)

计算结果如下所示：

> print(diff)
[1] 3.020833

现在我们开始计算p值，上面的代码中有一个diff变量，可以称为observed effect size，这是一个随机变量，我们还知道，在零假设的前提下，diff均值的分布是0，而这个随机变量分布的标准误则是总体的标准差(population standard deviation)除了样本数目的平方根，如下所示：
S E(\overline{X})=\sigma / \sqrt{N}
我们使用样本样本标准差(sample standard deviation)来对总体的标准差(population standard deviation)进行估计，在R中，也可以使用的sd()函数来计算标准误(SE)，如下所示：

sd(control)/sqrt(length(control))

结果如下所示：

> sd(control)/sqrt(length(control))
[1] 0.8725323

这个值就是样本均值(sample average)的SE，不过实际上，我们相想要的是diff的SE。我们前面已经知道了，统计学理论告诉我们，两个随机变量差的方差是原两个随机变量方差的和，因此我们计算一下方差以及平方根：

se <- sqrt(
var(treatment)/length(treatment) +
var(control)/length(control)
)
se

结果如下所示：

> se
[1] 1.469867

统计学理论还告诉我们，如果我们将一个随机变量除以其SE，我们就会得到一个SE为1的新的随机变量，如下所示：

tstat <- diff/se
tstat

结果如下所示：

> tstat <- diff/se
> tstat
[1] 2.055174

上面计算出来的tstat就是我们称的t统计量(t-statistic)。这是两个随机变量的比值，因此它也是一个随机变量。一旦我们知道了这个随机的分布，我们就很容易计算出其p值。

前面我们提到，CLT可以告诉我们，如果有大样本(large sample size)，那么样本的均值mean(treatment)和mean(control)都服从正态分布。统计学理论告诉我们，两个正态分布的随机变量的差值还服从正态分布，因此CLT也告诉我们，tstat接近于均值为0(零假设)，SD为1（我们除以它的SE）的正态分布。

现在我们要计算一下p值，此时我们需要问的一个问题就是：一个服从正态分布的随机变量有多大的概率会大于diff？R语言中的有一个pnorm()函数有解决这个问题。pnorm(a)会以返回一个值，这个值就是概率，它表示在一个标准正态分布中，一个随机变量低于a 的概率，具体这个函数的原理与用法，可以参考这篇笔记《R语言中dnorm, pnorm, qnorm与rnorm以及随机数》。

如果要计算大于a的概率，那么可以使用1-pnorm(a)，如果我们想知道看到一些像diff这种极端事件的概率：例如要么小于（更多的时候指的是负值）-abs(diff)，要么是大于abs(diff)，那么我们也可以过计算这两个尾部区域的概率：

righttail <- 1-pnorm(abs(tstat))
lefttail <- pnorm(-abs(tstat))
pval <- lefttail + righttail
print(pval)

计算结果如下所示：

> print(pval)
[1] 0.0398622

在这个案例中，p值小于0.05，根据传统的p值阈值0.05，我们就可以下给结论，这种差值有统计学意义。

现在我们就面临一个问题，CLT只有在大样本的前提下才有效，12个样本这个数目够么？根据经验法则，通常大于30个样本（这仅仅是经验法则），就能满足CLT。我们计算的这个p值，只有在假设成立的前提下才是有效的近似值，而实际情况并非如此，还好，除了CLT之外，我们还有另外一种选择。

t分布的实际计算

如果一个总体的分布是正态分布，那么我们在不需要CLT思想的前提下，就可以计算出t-statistic的精确分布。但是，如果我们只有小样本，我们就很难知道总体是否是正态的。但是，对于一些常见的统计量，例如体重，我们就可以猜测它的总体分布非常接近于正态分布，因此我们可以使用这种近似。此外，我们可以使用对样本画qq图，它会显示出样本的分布，如下所示：

library(rafalib)
mypar(1,2)
qqnorm(treatment)
qqline(treatment,col=2)
qqnorm(control)
qqline(control,col=2)
image

如果我们使用这种近似计算，那么从统计学理论角度来讲我们随机变量tstat的分布就服从一个t分布(-t-distribution)，这种分布比正态分布更复杂。t分布没有像正态分布那样的位置参数( location parameter )，正态分布的参数是自由度(degrees of freedom)，R中的t.test()函数可以用于计算相应的p值，如下所示：

t.test(treatment, control)

计算结果如下所示：

> t.test(treatment, control)

    Welch Two Sample t-test

data:  treatment and control
t = 2.0552, df = 20.236, p-value = 0.053
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.04296563  6.08463229
sample estimates:
mean of x mean of y 
 26.83417  23.81333

如果只想看p值，那么使用以下代码即可，如下所示：

result <- t.test(treatment,control)
result$p.value

结果如下：

> result <- t.test(treatment,control)
> result$p.value
[1] 0.05299888

从计算结果来看，这个p值有点大，这是因为在CLT近似的计算中，tstat的分布实际上是固定的（对于大样本来说，分母实际上就是固定的），而t分布的近似计算则是考虑到了分母（差值的标准误）是一个随机变量，是不固定的，样本数目越小，那分母的变化就越大。

这样来看的，我们使用了两种方法，得到了2个不同的p值，这在数据分析过程中很常见，因为我们是使用了不同的前提，不同的近似，因此会得到不同的结果。在后面讲到的功效计算(power calculation)中，我们会讲到I类错误，II类错误。在这里我们只说一下，使用CLT近似会错误地拒绝( incorrectly reject )零假设(假阳性)，而t分布则会错误地接受(inccorrectly acept)零假设(假阴性)。

在R中计算t检验如下所示：

library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
control <- filter(dat,Diet=="chow") %>% select(Bodyweight)
treatment <- filter(dat,Diet=="hf") %>% select(Bodyweight)
t.test(treatment,control)

结果如下所示：

> t.test(treatment,control)

    Welch Two Sample t-test

data:  treatment and control
t = 7.1932, df = 735.02, p-value = 1.563e-12
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 2.231533 3.906857
sample estimates:
mean of x mean of y 
 30.48201  27.41281

置信区间(Confidence Intervals)

在生命科学研究中，仅在研究结论中报道p值是不够的，因为统计学上的显著性无法保证科学研究上的意义。例如在对体重的差异进行统计时，使用大样本的情况下，即使是1毫克的差异，也有可能有统计学意义，不过这种1毫克的差异是一个重要的发现么？我们能否说明，某个因素导致了不到1%的体重变化？仅报道p值是无法提供一些有价值的信息的，也不会提供有关效应量(the effect size)的信息，效应量(effect size)就是我们前面提到的观察到的差异(observed difference)。有的时候，效应量会用除以对照组的均值，因此会表示为百分比。

为了纠正仅报道p值偏差，另外一个统计量就提了出来，即置信区间(confidence intervals)。置信区间包括了估计的效应量(estimated effect size)，以及与这个估计有关的不确定性(uncertainty)，现在我们使用这批小鼠数据来计算一下置信区间。

总体均值的置信区间

在我们展示如何过计算两组差值的置信区间之前，我们先来看一下如何计算对照组雌性小鼠的总体均值的置信区间。随后，我们会计算样本的置留区间，如下所示：

dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
chowPopulation <- dat[dat$Sex=="F" & dat$Diet=="chow",3]
mu_chow <- mean(chowPopulation)
print(mu_chow)

计算结果如下所示：

> print(mu_chow)
[1] 23.89338

现在计算的结果就是 总体均值\mu_{X}。

现在我们对这个结果进行估计，在实际分析过程中，我们并不会获取所有的总体，因此我们也像前面计算p值那样，我们来演示一下，通过抽样的方法来计算置信区间，现在我们先抽30个样本，如下所示：

N <- 30
chow <- sample(chowPopulation, N)
print(mean(chow))

计算结果如下所示：

> print(mean(chow))
[1] 23.95033

我们知道这个均值是一个随机变量，因此样本均值并不是一个非常好的估计。实际上，我们使用这个案例只是为了计算这个样本的均值，每次计算并非完全一样，如下所示：

> N <- 30
> chow <- sample(chowPopulation, N)
> print(mean(chow))
[1] 23.48467
> 
> N <- 30
> chow <- sample(chowPopulation, N)
> print(mean(chow))
[1] 25.22967
> 
> N <- 30
> chow <- sample(chowPopulation, N)
> print(mean(chow))
[1] 23.732

如上所示，每次抽样计算的均值都不一样。现在我们来了解一下置信区间，置信区间一种报道你数据中样本均值的一种统计学表示方法，它能明确地显示出随机变量的变异程度。

现在我们的的样本是30，CLT告诉我们，\bar{X}或chow组的均值(mean)服从一个均值为\mu_{X}（这个是chowPopulaiton这个总体的均值，即mean(chowPopulation)，这个均值服从均值为\mu_{X}，标准误s_{X} / \sqrt{N}的标准分布，计算过程如下所示：

se <- sd(chow)/sqrt(N)
print(se)

计算结果如下所示：

> se <- sd(chow)/sqrt(N)
> print(se)
[1] 0.597656

定义区间

95%置信区间是一个随机区间，这个区间有95%概率落在我们估计值的参数中。这里需要注意的是，说95%的随机区间将会包括真值，并不等于说真值有95%的概率落在我们的区间中。为了构建置信区间，我们需要注意，CLT告诉我们，\sqrt{N}\left(\overline{X}-\mu_{X}\right) / s_{X}服从均值为0，SD为1的正态分布，这个概率也就是以下事件的概率：
-2 \leq \sqrt{N}\left(\overline{X}-\mu_{X}\right) / s_{X} \leq 2
在R中的计算为：

pnorm(2) - pnorm(-2)

计算结果如下所示：

> pnorm(2) - pnorm(-2)
[1] 0.9544997

这个值与95%很接近，也就是qnorm(1-0.05/2)的这个值，现在我们把上面的公式变换一下，把\mu_{X}放中间，就成了如下的样子：
\overline{X}-2 s_{X} / \sqrt{N} \leq \mu_{X} \leq \overline{X}+2 s_{X} / \sqrt{N}

也就是上面这个值的概率是95%。需要注意的是，区间\overline{X} \pm 2 s_{X} / \sqrt{N}是两个边界，它们是随机的。此外，我们需要明确一下，置信区间的定义是，95%的随机区间（random intervals)包含真正的固定值\mu_{X}。对于一个计算结果来说，计算出来的这个特定的区间，它要么包括固定的总体均值\mu，要么不包括。

现在我们来模拟一下其中的思路，如下所示：

N <- 30
chow <- sample(chowPopulation, N)
se <- sd(chow)/sqrt(N)
Q <- qnorm(1 - 0.05/2)
interval <- c(mean(chow)-Q*se, mean(chow) + Q*se)
interval
interval[1] < mu_chow & interval[2] > mu_chow

结果如下所示：

> interval
[1] 22.17089 24.58044
> interval[1] < mu_chow & interval[2] > mu_chow
[1] TRUE

从上面我们可以看到，这个区间覆盖了总体均值\mu_{X}或mean(chowPopulation)。但是，如果我们再取另外一批样本的话，有可能不会覆盖总体均值\mu_{X}，但是，统计学理论告诉我们，经过均数次抽样后，这个样本覆盖均值的概率是95%。由于我们能够获取总体数据，现在我们来看一些样本：

library(rafalib)
B <- 250
mypar()
plot(mean(chowPopulation)+c(-7, 7), c(1,1), type="n",
     xlab="weight", ylab = "interval", ylim=c(1,B))
abline(v=mean(chowPopulation))
for (i in 1:B){
  chow <- sample(chowPopulation, N)
  sd <- sd(chow)/sqrt(N)
  interval <- c(mean(chow)-Q*se, mean(chow)+Q*se)
  covered <- 
    mean(chowPopulation) <= interval[2] & mean(chowPopulation) >= interval[1]
  color <- ifelse(covered, 1,2)
  lines(interval, c(i, i), col= color)
}
image

在这个案例中，我们模拟了250次随机抽样，计算了250次置信区间，其中颜色表示这个区间是否包含总体均值，其中绿色是包括的，棕色是不包括的。

上面这个案例的置信区间比较大（这个案例除以的是是\sqrt{5}，而非\sqrt{30}），我们可以看到很多区间并未覆盖\mu_{X}，这是因为CLT错误地告诉了我们(chow)的均值是近似地接近于正态分布，实际上这个数据是不太符合正态分布的（在接近\pm \infty的地方，它有一个比较扁平的尾巴）。这种错误会影响我们计算Q值，而这个Q值在计算的时，是以正态分布为前提,，使用qnorm()来计算的。在这种情况下，数据会更符合t分布，现在我们使用qt()函数来计算Q值，如下所示：

mypar()
plot(mean(chowPopulation)+ c(-7,7), c(1,1),type="n",
     xlab="weight",ylab="interval",ylim=c(1,B))
abline(v=mean(chowPopulation))
Q <- qt(1-0.05/2, df=4)
N <- 5
for (i in 1:B){
  chow <- sample(chowPopulation, N)
  se <- sd(chow)/sqrt(N)
  interval <- c(mean(chow)-Q*se, mean(chow)+Q*se)
  covered <- mean(chowPopulation) <= interval[2] & mean(chowPopulation) >= interval[1]
  color <- ifelse(covered, 1,2)
  lines(interval, c(i,i), col=color)
}
image

上面这个图中，我们使用了小样本来模拟了250次计算，得到了它的95%置信区间，这个置信区间是基于t分布进行计算的。与前面的那个案例相比，这个案例中置信区间更大，这是因为t分布的尾部更为扁平，得到的置信区间就更大，现在我们来计算一下t分布的累积分分布，并拿它与正态分布比较一下，如下所示：

qt(1-0.05/2, df=4)
qnorm(1-0.05/2)

计算结果如下所示：

> qt(1-0.05/2, df=4)
[1] 2.776445
> qnorm(1-0.05/2)
[1] 1.959964

从上面结果就可以看出来，t分布的置信区间更大，因此它才能以更大的概率(例如95%)来覆盖\mu_{X}。

置信区间与p值的关系

作者推荐报道实验结果时使用置信区间，而非仅用p值，对于某些原因，我们可能只需要提供p值即可，或者是提供0.05或0.01水平上的统计结果，其实置信区间也能提供这些信息。

当我们谈论t检验(t-test)的p值时，我们实际上是在谈论我们观察到的这两个样本的差值(\overline{Y}-\overline{X})以及更极端情况下差值的概率，就像是谈论两个总体之间的总体均值之差等于0这种极端情况下的概率。因此我们可以构建一个含有这个观察到的差值的置信区间，不过我们不再写为\overline{Y}-\overline{X}这种形成，而是写为一种新的形式，即d \equiv \overline{Y}-\overline{X}。

当我们来写一个差值的95%置信区间时，使用CLT则会写为d \pm 2 s_{d} / \sqrt{N}这种形式，这个结果并不包含0（假阳性）。因为置信区间不含0，这就暗示了，d-2 s_{d} / \sqrt{N}>0或d+2 s_{d} / \sqrt{N}<0，这种表达形式还说明，\sqrt{N} d / s_{d}>2或\sqrt{N} d / s_{d}<2。这说明，t统计量比2更极端，反过来又说明，p值必然小于0.05（为了更加精确地计算，可以使用qnorm(0.05/2)来替换2）。如果使用t分布来取代CLT，则使用使用qt(0.05/2, df=N-2。总之，一个95%或99%置信区间并不包含0，p它们的p值必然小于0.05或小于0.01。

现在我们使用t检验来计算一下d的置信区间，如下所示：

t.test(treatment,control,conf.level=0.95)$conf.int

结果如下所示：

> t.test(treatment,control,conf.level=0.95)$conf.int
[1] 2.231533 3.906857
attr(,"conf.level")
[1] 0.95

如果我们把置信水平改为0.9，如下所示：

> t.test(treatment,control,conf.level=0.90)$conf.int
[1] 2.366479 3.771911
attr(,"conf.level")
[1] 0.9

可以发现，置信区间变小了。

功效计算

在前面的案例中，我们研究了不同包含对小鼠体重的影响，由于我们在使用这个案例的时候，我们研究的是总体，发现了这两种包含对小鼠的体重有着明显的影响，计算结果如下所示：

library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% dplyr::select(Bodyweight) %>% unlist
hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% dplyr::select(Bodyweight) %>% unlist

mu_hf <- mean(hfPopulation)
mu_control <- mean(controlPopulation)
print(mu_hf - mu_control)
print((mu_hf - mu_control)/mu_control*100)

计算结果如下所示：

> print(mu_hf - mu_control)
[1] 2.375517
> print((mu_hf - mu_control)/mu_control*100)
[1] 9.942157

某些情况下，我们会采用抽样的形式，抽取一定数目的小鼠体重作为样本，此时就不能使用总体的计算方法，而是采用t检验，此时计算得到的p值不一定小于0.05，例如。我们随机抽取5只小鼠来计算一下，如下所示：

set.seed(1)
N <- 5
hf <- sample(hfPopulation, N)
control <- sample(controlPopulation, N)
t.test(hf,control)$p.value

计算结果如下所示：

> t.test(hf,control)$p.value
[1] 0.1410204

前后两种的计算方法出现了差异，第一次计算（使用小鼠的总体来进行计算）与第二次计算（随机抽两组的5只小鼠来计算）的p值不一样，是否出现了错误？

在第二次计算中，p值大于0.05，因此我们不拒绝零假设，从而说饮食对小鼠的体重没有影响么？能否这么讲？答案不是能。我们只能说，无法拒绝零假设（注意，没有后面那半句：包含对小鼠的体重没有影响）。我们无法拒绝零假设，并不意味着零假设一定为真。在这个特殊的案例中，这里我们面临的一个问题就是，我们没有足够的功效(power)。如果你从事科学研究，很有可能在某些时候不得不做一个功效计算。在很多情况下，这是一种首选义务，因此它可以帮助你在人的研究中，避免使用过多的小鼠，或者是避免使用过多的受试患者。现在我们就来解释一下统计功效(statistical power)的含义。

错误类型

我们要有这么一个意识，即，一旦我们做统计，就有可能犯错误，这也就是为什么p值不等于0。在进行统计时，当p值很小时，这说明我们可以拒绝零假设，但是这种拒绝也有可能犯错误（虽然犯错误的概率很低），例如零假设确实是真时。例如当p值等于0.05时，那么我们重复20实验，就有可能1次会发生零假设成功的事件，统计学家把这种错误称为I类错误(type I error)。

I类错误的定义就是，零假设成立时，我们拒绝了（其实不应该拒绝的），也就是所谓假阳性。这里为什么使用0.05，而不使用0.00001呢？这是因为如果阈值设得过低，成本会非常高。此时我们还要引入II类错误(type I error)，II类错误是指，零假设不成立，但是我们接受了零假设，也就是所谓的假阴性。前面我们抽样了5只小鼠，计算得到的值大于0.05，因此我们无法拒绝零假设（在0.05这个水平上），此时所犯的错误就是II类错误，也就是假阴性（从总体上看，实际上是有差异的，但是我们没有发现）。如果我们把阈值提高到0.25，我们就不会犯这个错误了（计算出的p值为0.1410204），不过通常情况下，不会这么干。

0.05与0.01阈值(cut-off)

多数的期刊或监管机构都坚持0.01或0.05的显著性水平，当然这么做无可厚非。我们本书的目标就之一就是让读者对p值以及置信区间有一个更深入的理解，以便读者能够在一些情况下做出合理的选择。不幸的是，这些阈值在某些程度上已经滥用，不过这是又是一个容易引起争论的话量，此处不表。

功效计算

功效(power)是指：当零假设为假时，拒绝它的概率，原书的描述如下所示：

Power is the probability of rejecting the null when the null is false.

不过“零假设为假”这种情况比较复杂，因为在许多情况下，零假设都有可能为假。\Delta \equiv \overline{Y}-\overline{X}可以是任意的，功效实际上要依赖于这个参数，功效同时还依赖于你估计值的标准误，反过来，标准误又依赖于样本数目，以及总体标准差。实际情况中，我们并不知道一些\Delta，\sigma_{X}，\sigma_{Y}值以及样本数目这些值。但是通过统计学理论，我们可以计算功效，R中的pwr包可以进行这样的计算，现在我们来模拟一下这个过程：

第一步：确定样本数目，我们确定为12，即N <- 12；

第二步：设定拒绝零假设的阈值，这里为0.05，也就是常见的p值水平，即alpha <- 0.05；

第三步：设计重复抽样次数，我们会重复地进行抽样，来计算一下拒绝零假设的比例，这里设为2000，即B <- 2000；

这个模拟过程就是：我们会从对照组和处理组中抽取N个样本，使用t检验来计算p值，观察p值是否小于0.05，现在我们把这个过程写为一个函数，代码如下所示：

N <- 12
alpha <- 0.05
B <- 2000

reject <- function(N, alpha=0.05){
  hf <- sample(hfPopulation, N)
  control <- sample(controlPopulation, N)
  pval <- t.test(hf, control)$p.value
  pval < alpha
}

如果样本数目为12，那么我们是否能拒绝零假设呢，如下所示：

> reject(12)
[1] FALSE

答案是不能。现在我们使用replicate()函数来重复2000次，如下所示：

N <- 12
rejections <- replicate(B, reject(N))
mean(rejections)

计算结果如下所示：

> mean(rejections)
[1] 0.218

其中rejecctions的结果是一个逻辑值，即里面只有TRUE或FALSE，其中TRUE为1，FALSE是0，那么我们计算的均值，即mean(rejections)就是TRUE所占的比例，因此我们计算出来的功效为0.218。这就是为什么，当我们知道零假设为假时，t检验也无法拒绝的原因。当样本数目为12时，功效是21.8%，对为能够在0.05水平上降低假阳性，那么我们可以设置更高的阈值，但这会导致II类错误的增多。

现在我们来看一下，功效是如何随着N的变化而变化的，如下所示：

Ns <- seq(5,50, 5)
power <- sapply(Ns, function(N){
  rejections <- replicate(B, reject(N))
  mean(rejections)
})
power

计算结果如下所示：

> power
 [1] 0.0960 0.1640 0.2775 0.3775 0.4795 0.5540 0.6415 0.7075 0.7725 0.8160

从上面的结果我们可以发现，随着N的增加，功效也在增加，曲线如下所示：

plot(Ns, power, type = "b")
image

再来看一个案例。我们把alpha这个拒绝阈值改变一下，再来看一下功效的变化。如果想要降低I类错误，那么功效就会越小，也就说，我们需要在两类错误之间进行权衡，现在我们看下面的代码，在保证N值相同的前提下，改变alpha后功效的大小：

N <- 30
alphas <- c(0.1, 0.05, 0.01, 0.001, 0.001)
power <- sapply(alphas, function(alpha){
  rejections <- replicate(B, reject(N, alpha=alpha))
  mean(rejections)
})
plot(alphas, power, xlab="alpha", type="b", log="x")

变化曲线如下所示：

image

在实际统计中，没有绝对的功效，也没有绝对的alpha值，但是我们要理解这背后的统计学硬。

替代假设下的p值

当零假设为假，替代假设为真时，那么p值有的时候就比较任意了。当替代假设为真时，我们通过增加样本数目，可以得到一个尽可能小的p值，通过从总体中不断抽取更大的样本数，我们就会知道p值的这种特性，现在看一下这个过程。

第一步：构建一个函数，这个函数能够返回一个样本数目为N时的p值，如下所示：

calculatePvalue <- function(N){
  hf <- sample(hfPopulation,N)
  control <- sample(controlPopulation, N)
  t.test(hf, control)$p.value
}

第二步：设定样本数目，对于每个样本数目的计算，我们会得到一系列p值，如下所示：

Ns <- seq(10, 200, by=10)
Ns_rep <- rep(Ns, each=10)

第三步：计算p值，如下所示：

pvalues <- sapply(Ns_rep, calculatePvalue)

第四步：绘制样本数目与其对应的p值，如下所示：

plot(Ns_rep, pvalues, log="y", xlab="Sample size",
     ylab="p-value")
abline(h=c(0.01, 0.05), col="red", lwd=2)
image

上图是随着样本数目的变化，p值的变化。从上面结果可以看出，当替代假设成立时，随着样本数目的增多，p值会越来越小，从上面的两条红线分别是0.01与0.05的区间。

一旦我们在某个阈值上拒绝了零假设，我们如果想要获得一个更小的p值，那么就需要更多的样本，例如实验小鼠。样本增大会增加我们对差值\Delta估计的精确性，但是，这个p值的降低仅仅是数学计算的自然结果。因为在计算p值时，它公式中的分母是样本数目的平方根，即\sqrt{N}，即如果\Delta不等于0，随着N的增加，p值必然降低。

因此，一个好的统计学结果需要列出效应量(effect size)与置信区间。效应量的计算是：将差值(difference)和置信区间(confidence interval)除以对照总体的均值，获得一个百分数，如下所示：

N <- 12
hf <- sample(hfPopulation, N)
control <- sample(controlPopulation, N)
diff <- mean(hf) - mean(control)
diff/mean(control)*100
t.test(hf, control)$conf.int / mean(control) * 100

计算结果为：

> diff/mean(control)*100
[1] 5.783149
> t.test(hf, control)$conf.int / mean(control) * 100
[1] -8.169363 19.735661
attr(,"conf.level")
[1] 0.95

此外，我们还可以写出一个名为Cohen's d的统计量，它是两组之间的差值除以这两组总和的标准差，如下所示：

sd_pool <- sqrt(((N-1)*var(hf) + (N-1)*var(control))/(2*N - 2))
diff / sd_pool

计算结果如下所示：

> sd_pool <- sqrt(((N-1)*var(hf) + (N-1)*var(control))/(2*N - 2))
> diff / sd_pool
[1] 0.3510624

这个结果告诉我们，hf组的小鼠体重的均值与对照组(control)小鼠体重平均值偏离多少个标准差。在替代假设下，t统计量会随着样本数目的增加而增加，但是效应量(effect size)与Cohen's d却会变得更加精确。

练习

P76

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 204,530评论 6赞 478
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 86,403评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 151,120评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,770评论 1赞 277
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,758评论 5赞 367
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,649评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,021评论 3赞 398
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,675评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,931评论 1赞 299
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,659评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,751评论 1赞 330
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,410评论 4赞 321
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,004评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,969评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1赞 260
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,042评论 2赞 350
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,493评论 2赞 343

推荐阅读更多精彩内容