title: DALS002-統(tǒng)計推斷(Inference)01-概率分布
date: 2019-07-22 12:0:00
type: "tags"
tags:

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

前言

這篇筆記是《Data Analysis for the Life Sciences》的第2章页滚，統(tǒng)計推斷(Inference)蛹含。這一部分的主要內(nèi)容涉及p值佳遂，置信區(qū)間(confidence intervals)台妆。

先來看一篇文獻，文獻信息如下所示：

Winzell, M. S. and B. Ahren (2004). "The high-fat diet-fed mouse: a model for studying mechanisms and treatment of impaired glucose tolerance and type 2 diabetes." Diabetes 53 Suppl 3: S215-219.
這篇文獻的摘要部分有這么一句話：

Body weight was higher in mice fed the high-fat diet already after the first week, due
to higher dietary intake in combination with lower metabolic efficiency.

上面摘要的這句話的意思是词身，飼喂高脂飲料的小鼠冶匹，在前一周內(nèi)，體重上升叙谨，這是由能量攝入較高，外加代謝消耗低導致的保屯。

為了支持這種結(jié)論手负，作者在結(jié)果部分提到了這些話：

Already during the first week after introduction of high-fat diet, body weight increased
significantly more in the high-fat diet-fed mice (+1.6± 0.1g) than in the normal diet-
fed mice (+ 0.2±0.1 g; P < 0.001)

當小鼠在第1周內(nèi)飼喂高脂飲料后，飲喂高脂飼料的小鼠體重的增加 (+1.6± 0.1g) 明顯高于普通飼料飼喂的小鼠體重增加(+0.2±0.1g; P < 0.001)姑尺。

這里面有幾個地方需要注意虫溜，第一就是P<0.001。第二就是±股缸。

要了解這兩個內(nèi)容衡楞，第一步就需要理解隨機變量(random variables)，這里使用作者在Github上提供的原始數(shù)據(jù)來計算敦姻，在第1篇筆記中已經(jīng)提到瘾境，導入數(shù)據(jù)如下所示：

filename <- file.path(dir,"extdata/femaleMiceWeights.csv")
dat <- read.csv(filename)

我們先要來研究一個問題：給小鼠特定的飼料，幾周后小鼠的體重是否會增加镰惦？現(xiàn)在我們從Jackson實驗室（Jackson實驗室是全球最大的實驗小鼠供應商）購買了24只小鼠迷守，隨機分布它們吃正常的飼料(chow)和高脂飼料(hf, high fat)，幾周后旺入，研究人員稱量小鼠的體重兑凿，獲取數(shù)據(jù)如下所示：

> head(dat)
  Diet Bodyweight
1 chow      21.51
2 chow      28.14
3 chow      24.04
4 chow      23.45
5 chow      23.68
6 chow      19.79

現(xiàn)在問題來了，飼喂高脂飼料(hf)的小鼠體重是否更重茵瘾？從數(shù)據(jù)中我們可以知道礼华，在hf組，編號24的小鼠體重是20.73g拗秘，它是最輕的小鼠圣絮，但是編號21的小鼠是34.02g，它是最重的小鼠雕旨，但是它們吃的都是高脂飼料扮匠，無法比較。此時凡涩，我們把這些小鼠體重的變化稱為變異(variability)棒搜。

再回到那篇文獻，那篇文獻下結(jié)論的時候活箕，指的變量就是小鼠體重的平均變化力麸，現(xiàn)在我們計算每組小鼠的平均變化，如下所示：

library(dplyr)
control <- filter(dat, Diet == "chow") %>% select(Bodyweight) %>% unlist
treatment <- filter(dat, Diet == "hf") %>% select(Bodyweight) %>% unlist
print(mean(treatment))
print(mean(control))
obsdiff <- mean(treatment) - mean(control)
obsdiff

結(jié)果如下所示：

> print(mean(treatment))
[1] 26.83417
> print(mean(control))
[1] 23.81333
> obsdiff <- mean(treatment) - mean(control)
> obsdiff
[1] 3.020833

從結(jié)果中我們可以發(fā)現(xiàn)，hf組小鼠的平均體重大概高于chow組小鼠平均體重的10%末盔。為什么在這個地方我們不能下確定的結(jié)論，還需要一個p值(p-value)與置信區(qū)間(confidence intervals)座慰？

這是因為這些均數(shù)數(shù)是隨機變量陨舱，如果我們再重復一次實驗，也就是說版仔，我們再從Jackson實驗室購買24只小鼠游盲，隨機將它們分為2組，一組hf蛮粮，一組chow益缎，我們還會得到一個不同的平均值。每重復一次然想，都會得到一個平均值莺奔，我們稱這樣的數(shù)據(jù)為隨機變量(random variable)。

隨機變量(random variable)

現(xiàn)在我們進一步研究隨機變量变泄。試想令哟，我們能夠獲取所有對照組小鼠的體重，把它們都放在R中計算妨蛹，從統(tǒng)計學角度來說屏富，這個所有的對照組小鼠的體重就是總體(population)，我們可以從這個總體中隨機抽取24只小鼠蛙卤，事實上來說狠半，我們不可能拿到總體，因為太大了颤难，不過我們可以模擬這個思路神年，現(xiàn)在下載這個數(shù)據(jù)，如下所示：

dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/femaleControlsPopulation.csv")
population <- read.csv(filename)
population <- unlist(population) # turn it into a numeric

現(xiàn)在我們隨機從總體(population)中抽取12個樣本行嗤，如下所示：

control <- sample(population, 12)
mean(control)
control <- sample(population, 12)
mean(control)
control <- sample(population, 12)
mean(control)

運行結(jié)果如下所示：

> control <- sample(population, 12)
> mean(control)
[1] 23.33917
> control <- sample(population, 12)
> mean(control)
[1] 23.53917
> control <- sample(population, 12)
> mean(control)
[1] 24.18167

現(xiàn)在我們可以發(fā)現(xiàn)瘤袖，每次抽取的時候，平均都發(fā)生了變化昂验，我們可以一直這么重復的抽樣下去捂敌，然后就能了解隨機變量的分布。

零假設(shè)(Null Hypothesis)

現(xiàn)在我們再回到前面的那個hf組與chow組小鼠平均體重差值的變量obsdiff上來既琴。眾所周知占婉，科學家都有懷疑精神，我們做了一次實驗甫恩，我們?nèi)绾沃?code>obsdiff這個是由飲食造成的呢逆济？如果我們給這2組小鼠都喂同樣的飼料，結(jié)果會是什么呢？萬一obsdiff這個值更大咋辦奖慌？統(tǒng)計學家們就把這種情況稱為零假設(shè)(Null Hypothesis)抛虫。這個零(null)指的就是懷疑的方向：我認為這2組小鼠的體重沒有差異的可能性(possibility)有多大(credence)。（注：原話是：we give credence to the possibility that there is no difference)简僧。

因為我們假設(shè)了我們能夠獲得總體建椰，我們實際上可能盡可能多地觀察到，當飼料沒對體重沒有影響時的差值〉郝恚現(xiàn)在我們假設(shè)隨機取了24只對照小鼠棉姐，把它們分為2組，也就是每組12只啦逆，現(xiàn)在我們觀察一下它們平均體重值的差值伞矩，如下所示：

control <- sample(population, 12) # control mice
treatment <- sample(population, 12)
# We assume high diet have no effect on weight
print(mean(treatment) - mean(control))

結(jié)果如下所示：

> print(mean(treatment) - mean(control))
[1] 0.2275

現(xiàn)在我們重復這個過程1萬次（最好不要使用for循環(huán)，這里可以使用replicate()函數(shù)）夏志，如下所示：

n <- 10000
null <- vector("numeric", n)
for (i in 1:n){
  control <- sample(population, 12)
  treatment <- sample(population, 12)
  null[i] <- mean(treatment) - mean(control)
}

其中變量null就是我們稱的零分布(null distribution)乃坤，它里面有10000個差值，現(xiàn)在我們來計算一下沟蔑，這里面的10000個差值與前面obsdiff的比較（原來過計算的obsdiff為3.020833）侥袜，看一下大于原來的3.020833這個數(shù)字的比例有多大，如下所示：

mean(null >= 3.020833)

如下所示：

> mean(null >= 3.020833)
[1] 0.013

這個比例是0.013溉贿，也就是1.3%枫吧，也就是說，假設(shè)飼料對小鼠的體重沒有影響宇色，我們重復1萬次實驗九杂，只有1.3%的比例，也就是說只有130次實驗能說明“飼料對小鼠的體重沒影響”這個結(jié)論宣蠕，另外9870次我們都得不出這個結(jié)論例隆，因此，即使作為懷疑性極強的科學們家抢蚀，也只能說镀层，當飼料對體重沒有影響時，我們只能看到1.3%的的可能性皿曲。這就是p值唱逢。

下面是通過replicate()函數(shù)來實現(xiàn)上述過程，如下所示：

null_f <- function(x){
  control <- sample(x, 12)
  treatment <- sample(x, 12)
  null <- mean(treatment) - mean(control)
  return(null)
}

null_result <- replicate(n, null_f(population))
mean(null_result>=3.020833)

結(jié)果如下所示：

> mean(null_result>=3.020833)
[1] 0.0123

1.23%屋休。

分布(distributions)

分布簡單來說就是坞古，對大量數(shù)據(jù)的精簡描述，例如你測量了一個總體中的所有男性的身高劫樟，現(xiàn)在你需要把這些數(shù)據(jù)向一個外星人描述痪枫，而這個外星人從來沒到過地球织堂，數(shù)據(jù)如下所示：

library(UsingR)
x <- father.son$fheight

注：這個包UsingR是《Using R for Introductory Statistics》 Second Edition的數(shù)據(jù)。其中father.son這個數(shù)據(jù)集含有的是父親的身高與兒子的身高（單位是英寸）奶陈，用于研究Pearson相關(guān)系數(shù)易阳。

現(xiàn)在向外星人描述比較簡單的方法就是把它們都列出來，現(xiàn)在我們隨機選擇10個數(shù)據(jù)吃粒，如下所示：

round(sample(x,10),10)

如下所示：

> round(sample(x,10),10)
 [1] 65.25267 62.41134 68.83427 66.69404 65.98544 69.60692 71.96564 68.11283 67.93353
[10] 67.41537

累積分布函數(shù)(Cumulative Distribution Function, CDF)

通過查看這些數(shù)字潦俺，我們了解一個大概，但是非常低效声搁。但是我們可以通過一個分布曲線(distribution)來畫一下這些數(shù)據(jù)，現(xiàn)在我們定義一下分布捕发，先取一個值a疏旨，小于a值的數(shù)字的比例，就是分布扎酷，如下所示：

這個函數(shù)叫累積分布函數(shù)(CDF)檐涝，當CDF來源于一個數(shù)據(jù)集中，而不是理論值時法挨，我們就稱它為經(jīng)驗累積分布函數(shù)(empirical CDF, ECDF)谁榜，我們可以畫出這個分布曲線，如下所示：

smallest <- floor(min(x))
largest <- ceiling(max(x))
values <- seq(smallest, largest, len=300)
heightecdf <- ecdf(x)
plot(values, heightecdf(values), type="l",
     xlab="a(Height in inches)", ylab="Pr(x <= a)")

image

直方圖(Histograms)

ecdf這個函數(shù)在R中并不常用凡纳，現(xiàn)在我們使用區(qū)間的形式來給出間內(nèi)的比例窃植，如下所示：

使用直方圖更加直觀，因為我們多數(shù)感興趣的是某個區(qū)間的高度荐糜，例如70到71這個區(qū)間的比例巷怜，而很少會去看低于70這個比例，直方圖如下所示：

hist(x)

image

可以對這個直方圖做進一步的修飾暴氏，如下所示：

bins <- seq(smallest, largest)
hist(x, breaks=bins, xlab="Height(in inches)", main="Adlt men heights")

image

這種圖比純列數(shù)字更加直觀延塑，從圖上我們大致就能判斷出高于72英寸的人大概有70個人。

概率分布(Probability Distribution)

分布的一個更加重要的功能是用于描述一個隨機變量的可能結(jié)果答渔。與一組固定的數(shù)據(jù)不同关带，我們實際上無法觀測到隨機變量的所有可能結(jié)果，因此沼撕，我們在這種情況下宋雏，不描述比例(proportions)，而是用于描述概率(probabilities)务豺。例如好芭，當我們從上面的那組身高數(shù)據(jù)中取機挑一個數(shù)字，它落在a與b區(qū)間內(nèi)的概率就可以用下面的公式表示：

上面公式中的X表示隨機變量冲呢，這個公式定義了隨機變量的概率分布舍败。例如，在這種情況下，如果我們知道了零假設(shè)下的小鼠體重均值的差值分布是真的邻薯，那么裙戏，我們就能在這個分布上找到與我們實驗中計算出的值對應值的概率，也就是p值厕诡。

在前面的案例中累榜，我們計算了10000個零假設(shè)下的隨機變量概率分布（這種重復抽樣的方法叫蒙特卡羅法,Monte Carlo），現(xiàn)在我們再重復一次灵嫌，不過這次有所不同壹罚，在每重復做一次實驗，我們就把這個點添加到圖片上寿羞，如下所示：

library(rafalib)
n <- 100
nullplot(-5, 5, 1, 30, xlab="Observed differences(grams)", ylab="Frequencey")
totals <- vector("numeric",11)
for (i in 1:n){
  control <- sample(population, 12)
  treatment <- sample(population, 12)
  nulldiff <- mean(treatment) - mean(control)
  j <- pmax(pmin(round(nulldiff)+6, 11),1)
  totals[j] <- totals[j]+1
  text(j-6, totals[j], pch=15, round(nulldiff,1),cex=1)
  if(i <15) Sys.sleep(1) # See values appear slowly
}

image

上面的圖其實就是一個直方圖猖凛，obsdiff最初的值是3.020833，從上面的圖我們可以知道绪穆，與這個值接近的數(shù)字很少辨泳。

從前面我們計算的null可以看出來，如下所示：

hist(null, freq=TRUE)
abline(v=obsdiff, col="red", lwd=2)

image

上面是我們手動計算Pr(a)的過程玖院，但多數(shù)情況下菠红，一些分布的Pr(a)公式都是已知的，例如正態(tài)分布难菌。

正態(tài)分布(Normal Distribution)

正態(tài)分布是最常見的分布之一试溯，normal就是正常，常見的意思郊酒，它的公式如下所示：

在R的pnorm()函數(shù)中，a通常是負無窮猎塞，b則是你要計算的參數(shù)试读，另外，要計算區(qū)間的概率荠耽，還需要兩個參數(shù)钩骇，分別是均值與標準差，即pnomr(x, mu, sigma)铝量，這個公式計算的就是小于x的這值的概率倘屹，現(xiàn)在回到前面的obsdiff(3.020833)這個變量上來，我們計算一下：

1 - pnorm(obsdiff, mean(null),sd(null))

如下所示：

> 1 - pnorm(obsdiff, mean(null),sd(null))
[1] 0.01377017

小結(jié)

計算hf組與chow組小鼠的體重差異的p值是非常簡單的慢叨。為了做這個計算纽匙，我們可以買來所有的小鼠，來重復的計算拍谐，找到零假設(shè)烛缔，但是不現(xiàn)實馏段，因此我們要使用統(tǒng)計推斷(Statistical Inference)的思想，這種思想可以允許我們使用少量的樣本践瓷，例如24只小鼠就能得到一個相對靠譜的結(jié)論院喜。

練習

原書P37