999热热精品店,成人精品在线99视频,久久久久超碰碰

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

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

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

中心極限定理(CLT)

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

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

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

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

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

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

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

也會(huì)得到 $Y-X$ 的均值，其實(shí)這也相當(dāng)于 $Y-X=Y+aX$ ，其中 $a=-1$ ，這其實(shí)也就是說 $Y-X$ 的均值為 $\mu_{Y}-\mu_{X}$ 。

如果說 $X$ 和 $Y$ 是相互獨(dú)立的樣本，那么 $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$ 這個(gè)差值的方差也是兩個(gè)樣本的方差和，這有點(diǎn)反直覺，我們需要知道的就是，X與Y這兩個(gè)樣本是相互獨(dú)立的，它們確實(shí)不相互影響。

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

通過這種近似計(jì)算，我們就能夠?qū)⒂?jì)算p值的過程簡(jiǎn)單化，因?yàn)槲覀冎廊魏沃翟跇?biāo)準(zhǔn)正態(tài)分布下的比例。例如，這些值中只有5%的比例會(huì)大于2（絕對(duì)值），計(jì)算過程如下所示：

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

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

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

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

但是，這還沒完。因?yàn)槲覀儾⒉恢揽傮w的標(biāo)準(zhǔn)差，也就是不知道 $\sigma_{X}$ 與 $\sigma_{Y}$ 。它們是未知的總體參數(shù)，但是，我們可以通過樣本的標(biāo)準(zhǔn)差對(duì)它們進(jìn)行估算，估算的總體標(biāo)準(zhǔn)差我們不能希臘字母 $\mu$ 來表示，只能用 $s$ 來表示，那么它們的估計(jì)值如下所示：
$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é)方面計(jì)算的原因，這里不表，記住就行。但是，為了更加直觀地說明一些東西，我們還是要涉及一點(diǎn)。試想，如果你有2個(gè)數(shù)，把它們畫在數(shù)軸上，它們之間的距離如果是L，那么中間的那一點(diǎn)就是它們的平均值，這個(gè)平均值離它們的距離分別都是0.5L。因此，從這個(gè)角度來說，你從這2個(gè)數(shù)中的一個(gè)就能夠獲取這個(gè)信息（這涉及自由度的問題），而不需要2個(gè)數(shù)，這點(diǎn)不怎么重要，重要的一點(diǎn)是，我們是用了 $s_{X}$ 和 $s_{Y}$ 來估計(jì)總體的均值，即 $\sigma_{X}$ 和 $\sigma_{Y}$ ，因?yàn)榭傮w均值未知，只能靠樣本均值來估計(jì)。

因此我們就把前面的那個(gè)公式：
$\frac{\overline{Y}-\overline{X}}{\sqrt{\frac{\sigma_{X}^{2}}{M}+\frac{\sigma_{Y}^{2}}{N}}}$
定義為下面的這個(gè)樣子：
$\frac{\overline{Y}-\overline{X}}{\sqrt{\frac{s_{X}^{2}}{M}+\frac{s_{Y}^{2}}{N}}}$
如果 $M=N$ ，那么公式就可以化簡(jiǎn)為如下所示：
$\sqrt{N} \frac{\overline{Y}-\overline{X}}{\sqrt{s_{X}^{2}+s_{Y}^{2}}}$
CLT告訴我們，當(dāng)M與N足夠大的時(shí)候，上面的這個(gè)隨機(jī)變量就會(huì)服從均值為0，標(biāo)準(zhǔn)差為1的正態(tài)分布。

t分布

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

現(xiàn)在我們以小鼠的表型數(shù)據(jù)為例說明一下，我們創(chuàng)建2向量，一個(gè)當(dāng)作對(duì)照總體，一個(gè)當(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}$ 的分布，而不是隨機(jī)變量 $\bar{Y}$ 的分布?，F(xiàn)在我們來看一下這兩個(gè)總體的分布：

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

image

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

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

image

練習(xí)

P47

中心極限定理的實(shí)際運(yùn)用

現(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)在計(jì)算一下種群的標(biāo)準(zhǔn)差，這里不要使用R的sd()函數(shù)，因此這個(gè)函數(shù)是計(jì)算樣本的標(biāo)準(zhǔn)差的，它會(huì)除以(n-1)，當(dāng)我們想要對(duì)總體進(jìn)行估計(jì)時(shí)，按以下代碼計(jì)算：

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)

計(jì)算結(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é)計(jì)算上是準(zhǔn)確的，我們不使用sd()或var()函數(shù)，我們可以使用rafalib包中的popvar()與popsd()函數(shù)，如下所示：

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

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

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

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

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理論對(duì)這些數(shù)據(jù)的近似值的適應(yīng)情況。如果這些數(shù)據(jù)非常接近正態(tài)分布，那么這些數(shù)據(jù)點(diǎn)就會(huì)落在一個(gè)直線上（這個(gè)直接是正態(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é)果最合適，這是因?yàn)樗目傮w相對(duì)最接近正態(tài)分布，平均值也最接近于正態(tài)分布。在實(shí)際中，我們可以計(jì)算這樣一個(gè)比值：除以估計(jì)的標(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$ 時(shí)，CLT并不會(huì)提供一個(gè)有用的估計(jì)，當(dāng) $N=12$ 時(shí)，在較大值的點(diǎn)方面，稍微有點(diǎn)偏離。當(dāng) $N=25$ 或 $N=50$ 時(shí)，所有的點(diǎn)基本上都在直線上。

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

練習(xí)

P52

t檢驗(yàn)的實(shí)際計(jì)算

現(xiàn)在我們來看一下如何通過計(jì)算得到p值，先來導(dǎo)入數(shù)據(jù)，在這一步中，我們要確定哪些是treatment組，哪些是control組，并且計(jì)算出它們的均值差異，如下所示：

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)

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

> print(diff)
[1] 3.020833

現(xiàn)在我們開始計(jì)算p值，上面的代碼中有一個(gè)diff變量，可以稱為observed effect size，這是一個(gè)隨機(jī)變量，我們還知道，在零假設(shè)的前提下，diff均值的分布是0，而這個(gè)隨機(jī)變量分布的標(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)來對(duì)總體的標(biāo)準(zhǔn)差(population standard deviation)進(jìn)行估計(jì)，在R中，也可以使用的sd()函數(shù)來計(jì)算標(biāo)準(zhǔn)誤(SE)，如下所示：

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

結(jié)果如下所示：

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

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

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

結(jié)果如下所示：

> se
[1] 1.469867

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

tstat <- diff/se
tstat

結(jié)果如下所示：

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

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

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

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

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

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

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

> print(pval)
[1] 0.0398622

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

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

t分布的實(shí)際計(jì)算

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

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

image

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

t.test(treatment, control)

計(jì)算結(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

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

這樣來看的，我們使用了兩種方法，得到了2個(gè)不同的p值，這在數(shù)據(jù)分析過程中很常見，因?yàn)槲覀兪鞘褂昧瞬煌那疤幔煌慕?，因此?huì)得到不同的結(jié)果。在后面講到的功效計(jì)算(power calculation)中，我們會(huì)講到I類錯(cuò)誤，II類錯(cuò)誤。在這里我們只說一下，使用CLT近似會(huì)錯(cuò)誤地拒絕( incorrectly reject )零假設(shè)(假陽性)，而t分布則會(huì)錯(cuò)誤地接受(inccorrectly acept)零假設(shè)(假陰性)。

在R中計(jì)算t檢驗(yàn)如下所示：

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

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

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

在我們展示如何過計(jì)算兩組差值的置信區(qū)間之前，我們先來看一下如何計(jì)算對(duì)照組雌性小鼠的總體均值的置信區(qū)間。隨后，我們會(huì)計(jì)算樣本的置留區(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)

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

> print(mu_chow)
[1] 23.89338

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

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

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

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

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

我們知道這個(gè)均值是一個(gè)隨機(jī)變量，因此樣本均值并不是一個(gè)非常好的估計(jì)。實(shí)際上，我們使用這個(gè)案例只是為了計(jì)算這個(gè)樣本的均值，每次計(jì)算并非完全一樣，如下所示：

> 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

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

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

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

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

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

定義區(qū)間

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

pnorm(2) - pnorm(-2)

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

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

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

也就是上面這個(gè)值的概率是95%。需要注意的是，區(qū)間 $\overline{X} \pm 2 s_{X} / \sqrt{N}$ 是兩個(gè)邊界，它們是隨機(jī)的。此外，我們需要明確一下，置信區(qū)間的定義是，95%的隨機(jī)區(qū)間（random intervals)包含真正的固定值 $\mu_{X}$ 。對(duì)于一個(gè)計(jì)算結(jié)果來說，計(jì)算出來的這個(gè)特定的區(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

從上面我們可以看到，這個(gè)區(qū)間覆蓋了總體均值 $\mu_{X}$ 或mean(chowPopulation)。但是，如果我們?cè)偃×硗庖慌鷺颖镜脑挘锌赡懿粫?huì)覆蓋總體均值 $\mu_{X}$ ，但是，統(tǒng)計(jì)學(xué)理論告訴我們，經(jīng)過均數(shù)次抽樣后，這個(gè)樣本覆蓋均值的概率是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

在這個(gè)案例中，我們模擬了250次隨機(jī)抽樣，計(jì)算了250次置信區(qū)間，其中顏色表示這個(gè)區(qū)間是否包含總體均值，其中綠色是包括的，棕色是不包括的。

上面這個(gè)案例的置信區(qū)間比較大（這個(gè)案例除以的是是 $\sqrt{5}$ ，而非 $\sqrt{30}$ ），我們可以看到很多區(qū)間并未覆蓋 $\mu_{X}$ ，這是因?yàn)镃LT錯(cuò)誤地告訴了我們(chow)的均值是近似地接近于正態(tài)分布，實(shí)際上這個(gè)數(shù)據(jù)是不太符合正態(tài)分布的（在接近 $\pm \infty$ 的地方，它有一個(gè)比較扁平的尾巴）。這種錯(cuò)誤會(huì)影響我們計(jì)算Q值，而這個(gè)Q值在計(jì)算的時(shí)，是以正態(tài)分布為前提,，使用qnorm()來計(jì)算的。在這種情況下，數(shù)據(jù)會(huì)更符合t分布，現(xiàn)在我們使用qt()函數(shù)來計(jì)算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

上面這個(gè)圖中，我們使用了小樣本來模擬了250次計(jì)算，得到了它的95%置信區(qū)間，這個(gè)置信區(qū)間是基于t分布進(jìn)行計(jì)算的。與前面的那個(gè)案例相比，這個(gè)案例中置信區(qū)間更大，這是因?yàn)閠分布的尾部更為扁平，得到的置信區(qū)間就更大，現(xiàn)在我們來計(jì)算一下t分布的累積分分布，并拿它與正態(tài)分布比較一下，如下所示：

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

計(jì)算結(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)系

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

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

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

現(xiàn)在我們使用t檢驗(yàn)來計(jì)算一下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ū)間變小了。

功效計(jì)算

在前面的案例中，我們研究了不同包含對(duì)小鼠體重的影響，由于我們?cè)谑褂眠@個(gè)案例的時(shí)候，我們研究的是總體，發(fā)現(xiàn)了這兩種包含對(duì)小鼠的體重有著明顯的影響，計(jì)算結(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)

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

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

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

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

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

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

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

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

錯(cuò)誤類型

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

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

0.05與0.01閾值(cut-off)

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

功效計(jì)算

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

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

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

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

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

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

這個(gè)模擬過程就是：我們會(huì)從對(duì)照組和處理組中抽取N個(gè)樣本，使用t檢驗(yàn)來計(jì)算p值，觀察p值是否小于0.05，現(xiàn)在我們把這個(gè)過程寫為一個(gè)函數(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)

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

> mean(rejections)
[1] 0.218

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

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

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

計(jì)算結(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

再來看一個(gè)案例。我們把alpha這個(gè)拒絕閾值改變一下，再來看一下功效的變化。如果想要降低I類錯(cuò)誤，那么功效就會(huì)越小，也就說，我們需要在兩類錯(cuò)誤之間進(jìn)行權(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

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

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

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

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

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

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

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

第三步：計(jì)算p值，如下所示：

pvalues <- sapply(Ns_rep, calculatePvalue)

第四步：繪制樣本數(shù)目與其對(duì)應(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í)，隨著樣本數(shù)目的增多，p值會(huì)越來越小，從上面的兩條紅線分別是0.01與0.05的區(qū)間。

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

因此，一個(gè)好的統(tǒng)計(jì)學(xué)結(jié)果需要列出效應(yīng)量(effect size)與置信區(qū)間。效應(yīng)量的計(jì)算是：將差值(difference)和置信區(qū)間(confidence interval)除以對(duì)照總體的均值，獲得一個(gè)百分?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

計(jì)算結(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

此外，我們還可以寫出一個(gè)名為Cohen's d的統(tǒng)計(jì)量，它是兩組之間的差值除以這兩組總和的標(biāo)準(zhǔn)差，如下所示：

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

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

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

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

練習(xí)

P76

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

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

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

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

中心極限定理(CLT)

t分布

練習(xí)

中心極限定理的實(shí)際運(yùn)用

練習(xí)

t檢驗(yàn)的實(shí)際計(jì)算

t分布的實(shí)際計(jì)算

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

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

定義區(qū)間

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

功效計(jì)算

錯(cuò)誤類型

0.05與0.01閾值(cut-off)

功效計(jì)算

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

練習(xí)

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

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

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

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

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

中心極限定理(CLT)

t分布

練習(xí)

中心極限定理的實(shí)際運(yùn)用

練習(xí)

t檢驗(yàn)的實(shí)際計(jì)算

t分布的實(shí)際計(jì)算

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

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

定義區(qū)間

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

功效計(jì)算

錯(cuò)誤類型

0.05與0.01閾值(cut-off)

功效計(jì)算

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

練習(xí)

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

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

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av