DALS004-統(tǒng)計推斷(Inference)02-CLT與T-test

title: DALS003-統(tǒng)計推斷(Inference)02-CLT與T-test
date: 2019-07-22 12:0:00
type: "tags"
tags:

  • 統(tǒng)計推斷
    categories:
  • 生物統(tǒng)計

中心極限定理(Central Limit Theorem)與t分布

中心極限定理(CLT)

當(dāng)樣本的數(shù)目非常大的，一個隨機樣本的均值\bar{Y}就會服從正態(tài)分布，這個正態(tài)分布的均值是總體均值(\mu_{Y})，標(biāo)準(zhǔn)差是總體的標(biāo)準(zhǔn)差\sigma_{Y}除以樣本的數(shù)目N的平方根。

如果我們將一個隨機變量減去一個常數(shù)，那么這個新生成的隨機變量就會領(lǐng)銜這個常數(shù)，從數(shù)學(xué)角度來說，如果X是一個服從均值為\mu的分布，a是一個常數(shù)，那么X-a的均值的分布的均值就是\mu-a。這個也適應(yīng)用乘法與標(biāo)準(zhǔn)差(SD)。

如果X是一個隨機變量，它的分布的均值\mu和標(biāo)準(zhǔn)差(SD)為\sigma，而a是一個常數(shù)，那么aX這個新的隨機變量分布的均值是a\mu，標(biāo)準(zhǔn)差為|a|\sigma。為了說明這個問題，我們就假設(shè)，每只小鼠的體重都減去10，那么平均體重也應(yīng)該會減去10，類似的，如果我們將小鼠的體重單位由克(g, gram)換成毫克(mg, milligram)，也就相當(dāng)于把原來的數(shù)字都乘以了1000，那么數(shù)字的擴散程度也會變大。

這就說明了，如果我們獲得了N個樣本，那么這個公式\frac{\overline{Y}-\mu}{\sigma_{Y} / \sqrt{N}}就接近于均值為0，標(biāo)準(zhǔn)差為1的正態(tài)分布。

現(xiàn)在我們感興趣的是兩個樣本的均值的差值。這個數(shù)學(xué)公式就比較有用了。

如果我們有兩個隨機變量，即X和Y，它們的總體均值與總體均值差為\mu_{X}和\mu_{Y}，\sigma_{X}和\sigma_{Y}，然后我們就會得到這些結(jié)果：

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

也會得到Y-X的均值，其實這也相當(dāng)于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這兩個樣本是相互獨立的，它們確實不相互影響。

從數(shù)學(xué)角度來理解上面的內(nèi)容對我們的研究目（兩個樣本的均值，以及均值的差值）。因此這兩個樣本都服從正態(tài)分布，它們的差值也服從正態(tài)分布，方差是兩個樣本總體的方差之和。在零假設(shè)（也就是說這兩個樣本所代表的兩個總體的均值沒有差異）成立的前提下，樣本均值\bar{Y}-\bar{X}的差值的分布接近于均值為0（也就是說沒有差異），標(biāo)準(zhǔn)差為\sqrt{\sigma_{X}^{2}+\sigma_{Y}^{2}} / \sqrt{N}的正態(tài)分布，也就是說下面的這個公式：
\frac{\overline{Y}-\overline{X}}{\sqrt{\frac{\sigma_{X}^{2}}{M}+\frac{\sigma_{Y}^{2}}{N}}}
的值接近于標(biāo)準(zhǔn)正態(tài)分布（均值為0，標(biāo)準(zhǔn)差為1）。

通過這種近似計算，我們就能夠?qū)⒂嬎鉷值的過程簡單化，因為我們知道任何值在標(biāo)準(zhǔn)正態(tài)分布下的比例。例如，這些值中只有5%的比例會大于2（絕對值），計算過程如下所示：

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

計算結(jié)果如下所示：

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

從上面結(jié)果我們就可以知道，只有5%的可能性會大于2，我們前面計算的obstdiff的值是3.020833，因此12只小鼠足夠了，不需要買更多的小鼠。

但是，這還沒完。因為我們并不知道總體的標(biāo)準(zhǔn)差，也就是不知道\sigma_{X}與\sigma_{Y}。它們是未知的總體參數(shù)，但是，我們可以通過樣本的標(biāo)準(zhǔn)差對它們進行估算，估算的總體標(biāo)準(zhǔn)差我們不能希臘字母\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，這有數(shù)學(xué)方面計算的原因，這里不表，記住就行。但是，為了更加直觀地說明一些東西，我們還是要涉及一點。試想，如果你有2個數(shù)，把它們畫在數(shù)軸上，它們之間的距離如果是L，那么中間的那一點就是它們的平均值，這個平均值離它們的距離分別都是0.5L。因此，從這個角度來說，你從這2個數(shù)中的一個就能夠獲取這個信息（這涉及自由度的問題），而不需要2個數(shù)，這點不怎么重要，重要的一點是，我們是用了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告訴我們，當(dāng)M與N足夠大的時候，上面的這個隨機變量就會服從均值為0，標(biāo)準(zhǔn)差為1的正態(tài)分布。

t分布

CLT的計算依賴于大量的樣本，也就是我們指的漸進性結(jié)果。當(dāng)CLT不適用時，還需要另外一種不依賴于漸進性結(jié)果的選擇。當(dāng)原始總體來源于一個變量，也就是Y時，如果這個Y是服從均值為0的標(biāo)準(zhǔn)正態(tài)分布，那么我們就可以計算下面這個變量的分布：
\sqrt{N} \frac{\overline{Y}}{s_{Y}}
這是兩個隨機變量的比值，因此這個比值不一定服從正態(tài)分布。事實上，這個公式的分母會偶然變小，因此就相應(yīng)地會增加觀察到更大值的概率。Gosset最初發(fā)現(xiàn)的這種分布，由于他開始在論文中以Student的名義投的稿，因此就稱這種分布為Student‘s t-distribution，也就是t分布。

現(xiàn)在我們以小鼠的表型數(shù)據(jù)為例說明一下，我們創(chuàng)建2向量，一個當(dāng)作對照總體，一個當(dāng)作高脂總體，如下所示：

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

這里需要注意提，我們假設(shè)的是的y_{1},y_{2},\dots,y_{n}的分布，而不是隨機變量\bar{Y}的分布。現(xiàn)在我們來看一下這兩個總體的分布：

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

我們可以使用qq圖來確認(rèn)一下上面的分布是比較接近正態(tài)分布的，這里我們只作初步了解，更詳細(xì)的了解在后面。如果在qq圖上，點大致都藻在一條直線附近，它們就屬于正態(tài)分布，如下所示：

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

練習(xí)

P47

中心極限定理的實際運用

現(xiàn)在我們使用上面的數(shù)據(jù)來看一下，中心極限定理是如何逼近樣本的均值的：

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

結(jié)果如下所示：

> 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

現(xiàn)在計算一下種群的標(biāo)準(zhǔn)差，這里不要使用R的sd()函數(shù)，因此這個函數(shù)是計算樣本的標(biāo)準(zhǔn)差的，它會除以(n-1)，當(dāng)我們想要對總體進行估計時，按以下代碼計算：

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)

計算結(jié)果如下所示：

> 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

因此，為了在數(shù)學(xué)計算上是準(zhǔn)確的，我們不使用sd()或var()函數(shù)，我們可以使用rafalib包中的popvar()與popsd()函數(shù)，如下所示：

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

這里需要注意，在實際統(tǒng)計中，我們并不會計算總體參數(shù)，我們通常是由樣本來估計總體參數(shù)，現(xiàn)在我們從總體中抽樣，如下所示：

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

CLT告訴我們，對于更大的N來說，這些樣本都近拉于正態(tài)分布。作為一種經(jīng)驗，這個N通常是大于30，這僅是一個經(jīng)驗法則，具體來說，它還需要取決于總體的分布。這里我們實際上可以查看一下近似值，并對各種N值進行計算，代碼如下所示：

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

現(xiàn)在使用qq圖看一下CLT理論對這些數(shù)據(jù)的近似值的適應(yīng)情況。如果這些數(shù)據(jù)非常接近正態(tài)分布，那么這些數(shù)據(jù)點就會落在一個直線上（這個直接是正態(tài)分布的分位數(shù)(quantiles))。越偏離，說明數(shù)據(jù)越不符合正態(tài)分布，如下所示：

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

從上面結(jié)果可以看出來，第3組數(shù)據(jù)的擬合結(jié)果最合適，這是因為它的總體相對最接近正態(tài)分布，平均值也最接近于正態(tài)分布。在實際中，我們可以計算這樣一個比值：除以估計的標(biāo)準(zhǔn)誤差來判斷，如下所示：

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

從上面可以看出來，當(dāng)N=3時，CLT并不會提供一個有用的估計，當(dāng)N=12時，在較大值的點方面，稍微有點偏離。當(dāng)N=25或N=50時，所有的點基本上都在直線上。

也就是說，在這個案例中，只要N=12就能證明CLT，就像前言提到的那樣，在多數(shù)情況下，這種模擬并不會很好。我們在這里只是通過這種模擬手段來說明CLT背后的思想與局限。

練習(xí)

P52

t檢驗的實際計算

現(xiàn)在我們來看一下如何通過計算得到p值，先來導(dǎo)入數(shù)據(jù)，在這一步中，我們要確定哪些是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)

計算結(jié)果如下所示：

> print(diff)
[1] 3.020833

現(xiàn)在我們開始計算p值，上面的代碼中有一個diff變量，可以稱為observed effect size，這是一個隨機變量，我們還知道，在零假設(shè)的前提下，diff均值的分布是0，而這個隨機變量分布的標(biāo)準(zhǔn)誤則是總體的標(biāo)準(zhǔn)差(population standard deviation)除了樣本數(shù)目的平方根，如下所示：
S E(\overline{X})=\sigma / \sqrt{N}
我們使用樣本樣本標(biāo)準(zhǔn)差(sample standard deviation)來對總體的標(biāo)準(zhǔn)差(population standard deviation)進行估計，在R中，也可以使用的sd()函數(shù)來計算標(biāo)準(zhǔn)誤(SE)，如下所示：

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

結(jié)果如下所示：

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

這個值就是樣本均值(sample average)的SE，不過實際上，我們相想要的是diff的SE。我們前面已經(jīng)知道了，統(tǒng)計學(xué)理論告訴我們，兩個隨機變量差的方差是原兩個隨機變量方差的和，因此我們計算一下方差以及平方根：

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

結(jié)果如下所示：

> se
[1] 1.469867

統(tǒng)計學(xué)理論還告訴我們，如果我們將一個隨機變量除以其SE，我們就會得到一個SE為1的新的隨機變量，如下所示：

tstat <- diff/se
tstat

結(jié)果如下所示：

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

上面計算出來的tstat就是我們稱的t統(tǒng)計量(t-statistic)。這是兩個隨機變量的比值，因此它也是一個隨機變量。一旦我們知道了這個隨機的分布，我們就很容易計算出其p值。

前面我們提到，CLT可以告訴我們，如果有大樣本(large sample size)，那么樣本的均值mean(treatment)和mean(control)都服從正態(tài)分布。統(tǒng)計學(xué)理論告訴我們，兩個正態(tài)分布的隨機變量的差值還服從正態(tài)分布，因此CLT也告訴我們，tstat接近于均值為0(零假設(shè))，SD為1（我們除以它的SE）的正態(tài)分布。

現(xiàn)在我們要計算一下p值，此時我們需要問的一個問題就是：一個服從正態(tài)分布的隨機變量有多大的概率會大于diff？R語言中的有一個pnorm()函數(shù)有解決這個問題。pnorm(a)會以返回一個值，這個值就是概率，它表示在一個標(biāo)準(zhǔn)正態(tài)分布中，一個隨機變量低于a 的概率，具體這個函數(shù)的原理與用法，可以參考這篇筆記《R語言中dnorm, pnorm, qnorm與rnorm以及隨機數(shù)》。

如果要計算大于a的概率，那么可以使用1-pnorm(a)，如果我們想知道看到一些像diff這種極端事件的概率：例如要么小于（更多的時候指的是負(fù)值）-abs(diff)，要么是大于abs(diff)，那么我們也可以過計算這兩個尾部區(qū)域的概率：

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

計算結(jié)果如下所示：

> print(pval)
[1] 0.0398622

在這個案例中，p值小于0.05，根據(jù)傳統(tǒng)的p值閾值0.05，我們就可以下給結(jié)論，這種差值有統(tǒng)計學(xué)意義。

現(xiàn)在我們就面臨一個問題，CLT只有在大樣本的前提下才有效，12個樣本這個數(shù)目夠么？根據(jù)經(jīng)驗法則，通常大于30個樣本（這僅僅是經(jīng)驗法則），就能滿足CLT。我們計算的這個p值，只有在假設(shè)成立的前提下才是有效的近似值，而實際情況并非如此，還好，除了CLT之外，我們還有另外一種選擇。

t分布的實際計算

如果一個總體的分布是正態(tài)分布，那么我們在不需要CLT思想的前提下，就可以計算出t-statistic的精確分布。但是，如果我們只有小樣本，我們就很難知道總體是否是正態(tài)的。但是，對于一些常見的統(tǒng)計量，例如體重，我們就可以猜測它的總體分布非常接近于正態(tài)分布，因此我們可以使用這種近似。此外，我們可以使用對樣本畫qq圖，它會顯示出樣本的分布，如下所示：

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

如果我們使用這種近似計算，那么從統(tǒng)計學(xué)理論角度來講我們隨機變量tstat的分布就服從一個t分布(-t-distribution)，這種分布比正態(tài)分布更復(fù)雜。t分布沒有像正態(tài)分布那樣的位置參數(shù)( location parameter )，正態(tài)分布的參數(shù)是自由度(degrees of freedom)，R中的t.test()函數(shù)可以用于計算相應(yīng)的p值，如下所示：

t.test(treatment, control)

計算結(jié)果如下所示：

> 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

結(jié)果如下：

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

從計算結(jié)果來看，這個p值有點大，這是因為在CLT近似的計算中，tstat的分布實際上是固定的（對于大樣本來說，分母實際上就是固定的），而t分布的近似計算則是考慮到了分母（差值的標(biāo)準(zhǔn)誤）是一個隨機變量，是不固定的，樣本數(shù)目越小，那分母的變化就越大。

這樣來看的，我們使用了兩種方法，得到了2個不同的p值，這在數(shù)據(jù)分析過程中很常見，因為我們是使用了不同的前提，不同的近似，因此會得到不同的結(jié)果。在后面講到的功效計算(power calculation)中，我們會講到I類錯誤，II類錯誤。在這里我們只說一下，使用CLT近似會錯誤地拒絕( incorrectly reject )零假設(shè)(假陽性)，而t分布則會錯誤地接受(inccorrectly acept)零假設(shè)(假陰性)。

在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)

結(jié)果如下所示：

> 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

置信區(qū)間(Confidence Intervals)

在生命科學(xué)研究中，僅在研究結(jié)論中報道p值是不夠的，因為統(tǒng)計學(xué)上的顯著性無法保證科學(xué)研究上的意義。例如在對體重的差異進行統(tǒng)計時，使用大樣本的情況下，即使是1毫克的差異，也有可能有統(tǒng)計學(xué)意義，不過這種1毫克的差異是一個重要的發(fā)現(xiàn)么？我們能否說明，某個因素導(dǎo)致了不到1%的體重變化？僅報道p值是無法提供一些有價值的信息的，也不會提供有關(guān)效應(yīng)量(the effect size)的信息，效應(yīng)量(effect size)就是我們前面提到的觀察到的差異(observed difference)。有的時候，效應(yīng)量會用除以對照組的均值，因此會表示為百分比。

為了糾正僅報道p值偏差，另外一個統(tǒng)計量就提了出來，即置信區(qū)間(confidence intervals)。置信區(qū)間包括了估計的效應(yīng)量(estimated effect size)，以及與這個估計有關(guān)的不確定性(uncertainty)，現(xiàn)在我們使用這批小鼠數(shù)據(jù)來計算一下置信區(qū)間。

總體均值的置信區(qū)間

在我們展示如何過計算兩組差值的置信區(qū)間之前，我們先來看一下如何計算對照組雌性小鼠的總體均值的置信區(qū)間。隨后，我們會計算樣本的置留區(qū)間，如下所示：

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)

計算結(jié)果如下所示：

> print(mu_chow)
[1] 23.89338

現(xiàn)在計算的結(jié)果就是 總體均值\mu_{X}。

現(xiàn)在我們對這個結(jié)果進行估計，在實際分析過程中，我們并不會獲取所有的總體，因此我們也像前面計算p值那樣，我們來演示一下，通過抽樣的方法來計算置信區(qū)間，現(xiàn)在我們先抽30個樣本，如下所示：

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

計算結(jié)果如下所示：

> 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

如上所示，每次抽樣計算的均值都不一樣?，F(xiàn)在我們來了解一下置信區(qū)間，置信區(qū)間一種報道你數(shù)據(jù)中樣本均值的一種統(tǒng)計學(xué)表示方法，它能明確地顯示出隨機變量的變異程度。

現(xiàn)在我們的的樣本是30，CLT告訴我們，\bar{X}或chow組的均值(mean)服從一個均值為\mu_{X}（這個是chowPopulaiton這個總體的均值，即mean(chowPopulation)，這個均值服從均值為\mu_{X}，標(biāo)準(zhǔn)誤s_{X} / \sqrt{N}的標(biāo)準(zhǔn)分布，計算過程如下所示：

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

計算結(jié)果如下所示：

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

定義區(qū)間

95%置信區(qū)間是一個隨機區(qū)間，這個區(qū)間有95%概率落在我們估計值的參數(shù)中。這里需要注意的是，說95%的隨機區(qū)間將會包括真值，并不等于說真值有95%的概率落在我們的區(qū)間中。為了構(gòu)建置信區(qū)間，我們需要注意，CLT告訴我們，\sqrt{N}\left(\overline{X}-\mu_{X}\right) / s_{X}服從均值為0，SD為1的正態(tài)分布，這個概率也就是以下事件的概率：
-2 \leq \sqrt{N}\left(\overline{X}-\mu_{X}\right) / s_{X} \leq 2
在R中的計算為：

pnorm(2) - pnorm(-2)

計算結(jié)果如下所示：

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

這個值與95%很接近，也就是qnorm(1-0.05/2)的這個值，現(xiàn)在我們把上面的公式變換一下，把\mu_{X}放中間，就成了如下的樣子：
\overline{X}-2 s_{X} / \sqrt{N} \leq \mu_{X} \leq \overline{X}+2 s_{X} / \sqrt{N}

也就是上面這個值的概率是95%。需要注意的是，區(qū)間\overline{X} \pm 2 s_{X} / \sqrt{N}是兩個邊界，它們是隨機的。此外，我們需要明確一下，置信區(qū)間的定義是，95%的隨機區(qū)間（random intervals)包含真正的固定值\mu_{X}。對于一個計算結(jié)果來說，計算出來的這個特定的區(qū)間，它要么包括固定的總體均值\mu，要么不包括。

現(xiàn)在我們來模擬一下其中的思路，如下所示：

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

結(jié)果如下所示：

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

從上面我們可以看到，這個區(qū)間覆蓋了總體均值\mu_{X}或mean(chowPopulation)。但是，如果我們再取另外一批樣本的話，有可能不會覆蓋總體均值\mu_{X}，但是，統(tǒng)計學(xué)理論告訴我們，經(jīng)過均數(shù)次抽樣后，這個樣本覆蓋均值的概率是95%。由于我們能夠獲取總體數(shù)據(jù)，現(xiàn)在我們來看一些樣本：

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次置信區(qū)間，其中顏色表示這個區(qū)間是否包含總體均值，其中綠色是包括的，棕色是不包括的。

上面這個案例的置信區(qū)間比較大（這個案例除以的是是\sqrt{5}，而非\sqrt{30}），我們可以看到很多區(qū)間并未覆蓋\mu_{X}，這是因為CLT錯誤地告訴了我們(chow)的均值是近似地接近于正態(tài)分布，實際上這個數(shù)據(jù)是不太符合正態(tài)分布的（在接近\pm \infty的地方，它有一個比較扁平的尾巴）。這種錯誤會影響我們計算Q值，而這個Q值在計算的時，是以正態(tài)分布為前提,，使用qnorm()來計算的。在這種情況下，數(shù)據(jù)會更符合t分布，現(xiàn)在我們使用qt()函數(shù)來計算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%置信區(qū)間，這個置信區(qū)間是基于t分布進行計算的。與前面的那個案例相比，這個案例中置信區(qū)間更大，這是因為t分布的尾部更為扁平，得到的置信區(qū)間就更大，現(xiàn)在我們來計算一下t分布的累積分分布，并拿它與正態(tài)分布比較一下，如下所示：

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

計算結(jié)果如下所示：

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

從上面結(jié)果就可以看出來，t分布的置信區(qū)間更大，因此它才能以更大的概率(例如95%)來覆蓋\mu_{X}。

置信區(qū)間與p值的關(guān)系

作者推薦報道實驗結(jié)果時使用置信區(qū)間，而非僅用p值，對于某些原因，我們可能只需要提供p值即可，或者是提供0.05或0.01水平上的統(tǒng)計結(jié)果，其實置信區(qū)間也能提供這些信息。

當(dāng)我們談?wù)搕檢驗(t-test)的p值時，我們實際上是在談?wù)撐覀冇^察到的這兩個樣本的差值(\overline{Y}-\overline{X})以及更極端情況下差值的概率，就像是談?wù)搩蓚€總體之間的總體均值之差等于0這種極端情況下的概率。因此我們可以構(gòu)建一個含有這個觀察到的差值的置信區(qū)間，不過我們不再寫為\overline{Y}-\overline{X}這種形成，而是寫為一種新的形式，即d \equiv \overline{Y}-\overline{X}。

當(dāng)我們來寫一個差值的95%置信區(qū)間時，使用CLT則會寫為d \pm 2 s_u0z1t8os / \sqrt{N}這種形式，這個結(jié)果并不包含0（假陽性）。因為置信區(qū)間不含0，這就暗示了，d-2 s_u0z1t8os / \sqrt{N}>0或d+2 s_u0z1t8os / \sqrt{N}<0，這種表達形式還說明，\sqrt{N} d / s_u0z1t8os>2或\sqrt{N} d / s_u0z1t8os<2。這說明，t統(tǒng)計量比2更極端，反過來又說明，p值必然小于0.05（為了更加精確地計算，可以使用qnorm(0.05/2)來替換2）。如果使用t分布來取代CLT，則使用使用qt(0.05/2, df=N-2?？傊?，一個95%或99%置信區(qū)間并不包含0，p它們的p值必然小于0.05或小于0.01。

現(xiàn)在我們使用t檢驗來計算一下d的置信區(qū)間，如下所示：

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

結(jié)果如下所示：

> 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

可以發(fā)現(xiàn)，置信區(qū)間變小了。

功效計算

在前面的案例中，我們研究了不同包含對小鼠體重的影響，由于我們在使用這個案例的時候，我們研究的是總體，發(fā)現(xiàn)了這兩種包含對小鼠的體重有著明顯的影響，計算結(jié)果如下所示：

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)

計算結(jié)果如下所示：

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

某些情況下，我們會采用抽樣的形式，抽取一定數(shù)目的小鼠體重作為樣本，此時就不能使用總體的計算方法，而是采用t檢驗，此時計算得到的p值不一定小于0.05，例如。我們隨機抽取5只小鼠來計算一下，如下所示：

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

計算結(jié)果如下所示：

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

前后兩種的計算方法出現(xiàn)了差異，第一次計算（使用小鼠的總體來進行計算）與第二次計算（隨機抽兩組的5只小鼠來計算）的p值不一樣，是否出現(xiàn)了錯誤？

在第二次計算中，p值大于0.05，因此我們不拒絕零假設(shè)，從而說飲食對小鼠的體重沒有影響么？能否這么講？答案不是能。我們只能說，無法拒絕零假設(shè)（注意，沒有后面那半句：包含對小鼠的體重沒有影響）。我們無法拒絕零假設(shè)，并不意味著零假設(shè)一定為真。在這個特殊的案例中，這里我們面臨的一個問題就是，我們沒有足夠的功效(power)。如果你從事科學(xué)研究，很有可能在某些時候不得不做一個功效計算。在很多情況下，這是一種首選義務(wù)，因此它可以幫助你在人的研究中，避免使用過多的小鼠，或者是避免使用過多的受試患者?，F(xiàn)在我們就來解釋一下統(tǒng)計功效(statistical power)的含義。

錯誤類型

我們要有這么一個意識，即，一旦我們做統(tǒng)計，就有可能犯錯誤，這也就是為什么p值不等于0。在進行統(tǒng)計時，當(dāng)p值很小時，這說明我們可以拒絕零假設(shè)，但是這種拒絕也有可能犯錯誤（雖然犯錯誤的概率很低），例如零假設(shè)確實是真時。例如當(dāng)p值等于0.05時，那么我們重復(fù)20實驗，就有可能1次會發(fā)生零假設(shè)成功的事件，統(tǒng)計學(xué)家把這種錯誤稱為I類錯誤(type I error)。

I類錯誤的定義就是，零假設(shè)成立時，我們拒絕了（其實不應(yīng)該拒絕的），也就是所謂假陽性。這里為什么使用0.05，而不使用0.00001呢？這是因為如果閾值設(shè)得過低，成本會非常高。此時我們還要引入II類錯誤(type I error)，II類錯誤是指，零假設(shè)不成立，但是我們接受了零假設(shè)，也就是所謂的假陰性。前面我們抽樣了5只小鼠，計算得到的值大于0.05，因此我們無法拒絕零假設(shè)（在0.05這個水平上），此時所犯的錯誤就是II類錯誤，也就是假陰性（從總體上看，實際上是有差異的，但是我們沒有發(fā)現(xiàn)）。如果我們把閾值提高到0.25，我們就不會犯這個錯誤了（計算出的p值為0.1410204），不過通常情況下，不會這么干。

0.05與0.01閾值(cut-off)

多數(shù)的期刊或監(jiān)管機構(gòu)都堅持0.01或0.05的顯著性水平，當(dāng)然這么做無可厚非。我們本書的目標(biāo)就之一就是讓讀者對p值以及置信區(qū)間有一個更深入的理解，以便讀者能夠在一些情況下做出合理的選擇。不幸的是，這些閾值在某些程度上已經(jīng)濫用，不過這是又是一個容易引起爭論的話量，此處不表。

功效計算

功效(power)是指：當(dāng)零假設(shè)為假時，拒絕它的概率，原書的描述如下所示：

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

不過“零假設(shè)為假”這種情況比較復(fù)雜，因為在許多情況下，零假設(shè)都有可能為假。\Delta \equiv \overline{Y}-\overline{X}可以是任意的，功效實際上要依賴于這個參數(shù)，功效同時還依賴于你估計值的標(biāo)準(zhǔn)誤，反過來，標(biāo)準(zhǔn)誤又依賴于樣本數(shù)目，以及總體標(biāo)準(zhǔn)差。實際情況中，我們并不知道一些\Delta，\sigma_{X}，\sigma_{Y}值以及樣本數(shù)目這些值。但是通過統(tǒng)計學(xué)理論，我們可以計算功效，R中的pwr包可以進行這樣的計算，現(xiàn)在我們來模擬一下這個過程：

第一步：確定樣本數(shù)目，我們確定為12，即N <- 12；

第二步：設(shè)定拒絕零假設(shè)的閾值，這里為0.05，也就是常見的p值水平，即alpha <- 0.05；

第三步：設(shè)計重復(fù)抽樣次數(shù)，我們會重復(fù)地進行抽樣，來計算一下拒絕零假設(shè)的比例，這里設(shè)為2000，即B <- 2000；

這個模擬過程就是：我們會從對照組和處理組中抽取N個樣本，使用t檢驗來計算p值，觀察p值是否小于0.05，現(xiàn)在我們把這個過程寫為一個函數(shù)，代碼如下所示：

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
}

如果樣本數(shù)目為12，那么我們是否能拒絕零假設(shè)呢，如下所示：

> reject(12)
[1] FALSE

答案是不能?，F(xiàn)在我們使用replicate()函數(shù)來重復(fù)2000次，如下所示：

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

計算結(jié)果如下所示：

> mean(rejections)
[1] 0.218

其中rejecctions的結(jié)果是一個邏輯值，即里面只有TRUE或FALSE，其中TRUE為1，FALSE是0，那么我們計算的均值，即mean(rejections)就是TRUE所占的比例，因此我們計算出來的功效為0.218。這就是為什么，當(dāng)我們知道零假設(shè)為假時，t檢驗也無法拒絕的原因。當(dāng)樣本數(shù)目為12時，功效是21.8%，對為能夠在0.05水平上降低假陽性，那么我們可以設(shè)置更高的閾值，但這會導(dǎo)致II類錯誤的增多。

現(xiàn)在我們來看一下，功效是如何隨著N的變化而變化的，如下所示：

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

計算結(jié)果如下所示：

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

從上面的結(jié)果我們可以發(fā)現(xiàn)，隨著N的增加，功效也在增加，曲線如下所示：

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

再來看一個案例。我們把alpha這個拒絕閾值改變一下，再來看一下功效的變化。如果想要降低I類錯誤，那么功效就會越小，也就說，我們需要在兩類錯誤之間進行權(quán)衡，現(xiàn)在我們看下面的代碼，在保證N值相同的前提下，改變alpha后功效的大?。?/p> 

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

在實際統(tǒng)計中，沒有絕對的功效，也沒有絕對的alpha值，但是我們要理解這背后的統(tǒng)計學(xué)硬。

替代假設(shè)下的p值

當(dāng)零假設(shè)為假，替代假設(shè)為真時，那么p值有的時候就比較任意了。當(dāng)替代假設(shè)為真時，我們通過增加樣本數(shù)目，可以得到一個盡可能小的p值，通過從總體中不斷抽取更大的樣本數(shù)，我們就會知道p值的這種特性，現(xiàn)在看一下這個過程。

第一步：構(gòu)建一個函數(shù)，這個函數(shù)能夠返回一個樣本數(shù)目為N時的p值，如下所示：

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

第二步：設(shè)定樣本數(shù)目，對于每個樣本數(shù)目的計算，我們會得到一系列p值，如下所示：

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

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

pvalues <- sapply(Ns_rep, calculatePvalue)

第四步：繪制樣本數(shù)目與其對應(yīng)的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

上圖是隨著樣本數(shù)目的變化，p值的變化。從上面結(jié)果可以看出，當(dāng)替代假設(shè)成立時，隨著樣本數(shù)目的增多，p值會越來越小，從上面的兩條紅線分別是0.01與0.05的區(qū)間。

一旦我們在某個閾值上拒絕了零假設(shè)，我們?nèi)绻胍@得一個更小的p值，那么就需要更多的樣本，例如實驗小鼠。樣本增大會增加我們對差值\Delta估計的精確性，但是，這個p值的降低僅僅是數(shù)學(xué)計算的自然結(jié)果。因為在計算p值時，它公式中的分母是樣本數(shù)目的平方根，即\sqrt{N}，即如果\Delta不等于0，隨著N的增加，p值必然降低。

因此，一個好的統(tǒng)計學(xué)結(jié)果需要列出效應(yīng)量(effect size)與置信區(qū)間。效應(yīng)量的計算是：將差值(difference)和置信區(qū)間(confidence interval)除以對照總體的均值，獲得一個百分?jǐn)?shù)，如下所示：

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

計算結(jié)果為：

> 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的統(tǒng)計量，它是兩組之間的差值除以這兩組總和的標(biāo)準(zhǔn)差，如下所示：

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

計算結(jié)果如下所示：

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

這個結(jié)果告訴我們，hf組的小鼠體重的均值與對照組(control)小鼠體重平均值偏離多少個標(biāo)準(zhǔn)差。在替代假設(shè)下，t統(tǒng)計量會隨著樣本數(shù)目的增加而增加，但是效應(yīng)量(effect size)與Cohen's d卻會變得更加精確。

練習(xí)

P76

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容