1、隨機變量( random variable)概念的引入
該數據來自傑克遜實驗室。2組數據，每組12只老鼠，一組普通食物，另一組高脂肪(hf)飲食。幾周后，科學家們稱了每只老鼠的體重，得到了這個數據:
dir <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/"
filename <- "femaleMiceWeights.csv"
url <- paste0(dir, filename)
dat <- read.csv(url)
library(dplyr)
control <- filter(dat,Diet=="chow") %>% select(Bodyweight) %>% unlist #其中%>%相當於管道符，fileter相當於Excel中按關鍵詞行篩選，select為列篩選，只保留你提到的變量
treatment <- filter(dat,Diet=="hf") %>% select(Bodyweight) %>% unlist
print(mean(treatment))
print(mean(control) )
obsdiff <- mean(treatment) - mean(control) #3.020833
print(obsdiff)
因此，食用hf的小鼠體重增加了10%我們做了什么?為什么我們需要p值和置信區間?原因是這些平均值是隨機變量。它們可以取很多值。如果我們重復這個實驗，我們從傑克遜實驗室得到24只新老鼠，然后隨機分配給每一種食物，我們得到不同的平均值。每次我們重復這個實驗，我們得到不同的值。我們稱之為隨機變量。
2、理解隨機變量( random variable)
population <- read.csv('https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/femaleControlsPopulation.csv')
str(population)
population <- unlist(population)
control <- sample(population,12)
mean(control)
control <- sample(population,12)
mean(control)
control <- sample(population,12)
mean(control)
可以看到每次抽樣12只，平均體重是有差別的。
3、 The Null Hypothesis
現在讓我們回到obsdiff的平均差異。作為懷疑論者，我們怎么知道這種強迫症是由飲食引起的呢?如果我們給24只老鼠同樣的食物會發生什么?我們會看到這么大的不同嗎?統計學家將這種情況稱為零假設。“null”這個名字是用來提醒我們，我們在扮演懷疑論者的角色:我們相信沒有區別的可能性。因為我們有群體( population)，我們可以隨機從群體中抽取24只老鼠，給以同樣的食物，然后將24只老鼠隨機分配兩組，來觀察這兩組的平均值。
Because we have access to the population, we can actually observe as many values as we want of the difference of the averages when the diet has no effect. We can do this by randomly sampling 24 control mice, giving them the same diet, and then recording the difference in mean between two randomly split groups of 12 and 12. Here is this process written in R code:
control <- sample(population,12)
treatment <- sample(population,12)
print(mean(treatment) - mean(control)) #注意每次結果會不一樣，因為是隨機抽取的
現在我們用for循環隨機抽取10000次。
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) }
mean(null >= obsdiff) ##0.0138
The values in null form what we call the null distribution.So what percent of the 10,000 are bigger than obsdiff?
在10000次模擬中，只有很小的一部分。作為懷疑論者，我們的結論是什么?在沒有飲食影響的情況下，這里隨機抽樣10000次和我們在實驗中觀察到的差異（框1）只有1.5%,這就是所謂的p值，即這1.5%的差異是隨機(偶然)產生的，該值越小，說明實驗值與對照值差別越大，處理效果越明顯(p越小，越拒絕零 假設)；該值越大說明實驗值與對照值差異是偶然因素造成的，越接受零假設。我們將在本書后面更正式地定義它。
Only a small percent of the 10,000 simulations. As skeptics what do we conclude? When there is no diet effect, we see a difference as big as the one we observed only 1.5% of the time. This is what is known as a p-value
4、distribution
我們已經解釋了零假設中零的含義，但分布到底是什么呢?考慮分布最簡單的方法是用許多數字的緊湊描述。例如，假設你測量了一個總體中所有人的身高。想象一下，你需要向一個不知道這些高度的人描述這些數字，比如一個從未去過地球的外星人。假設所有這些高度都包含在下面的數據集中，總結這些數字的一種方法是簡單地把它們列出來讓外星人看:
data(father.son,package="UsingR")
x <- father.son$fheight
round(sample(x,10),1)
5、Cumulative Distribution Function（累積分布函數）
通過瀏覽這些數字，我們開始大致了解整個列表的樣子，但它肯定是低效的。我們可以通過定義和可視化分布來快速改進這種方法。為了定義一個分布，我們對a的所有可能值，計算列表中低於a的數字的比例。我們使用以下符號:
F(a) =Pr(x ≤ a)
這稱為累積分布函數(cumulative distribution function ，CDF)。當CDF來自於數據時，與理論相對，我們也稱它為經驗CDF (ECDF)。
6、Histograms
雖然實證的CDF概念在統計教科書中得到了廣泛的討論，但在實際操作中卻並不十分流行。原因是直方圖給了我們相同的信息，而且更容易解釋。直方圖告訴我們值在區間中的比例:
Pr(a ≤ x ≤ b) = F (b) − F (a)。把這些高度用柱狀圖表示就是我們所說的直方圖。這是一個更有用的圖，因為我們通常更感興趣的是區間，在70到71英寸之間的這個百分比，而不是小於特定高度的那個百分比。通過直方圖也更容易區分分布的不同類型(族)。這是高度的直方圖:
hist(x,xlab='height',ylab = 'frequency',main='histogram')
我們可以通過以下方式來指定區間(bin)和添加更好的標簽:
smallest=min(x)
largest=max(x)
bins <- seq(smallest, largest+1)
hist(x,breaks=bins,xlab="Height (in inches)",main="Adult men heights")
向外星人展示這個情節比展示數字更有信息。用這個簡單的圖，我們可以近似出任意給定區間內的個體數量。例如，大約有70個人身高超過6英尺(72英寸)。
7、Probability Distribution
另外一個重要的概念是概率分布，我們不是描述比例，而是描述概率。例如，如果我們從列表中隨機選擇一個高度，那么它落在a和b之間的概率表示為：
注意，現在將X大寫，以區分它是一個隨機變量，上面的方程定義了隨機變量的概率分布。知道這個分布在科學上非常有用。例如，在上面的例子中，如果我們知道零假設為真時小鼠體重均值差的分布，也就是零分布（null distribution），我們就可以計算出觀察到一個和我們之前實驗處理組得到的一樣大的值的概率，也就是p值。
Knowing this distribution is incredibly useful in science. For example, in the case above, if we know the distribution of the difference in mean of mouse weights when the null hypothesis is true,referred to as the null distribution, we can compute the probability of observing a value as large as we did, referred to as a p-value.
讓我們重復上面的循環，但這次讓我們在每次重新運行實驗時向圖中添加一個點。在上一部分中，我們運行了所謂的MonteCarlo模擬(我們將在后面的章節中提供更多關於蒙特卡羅模擬的細節)，在零假設下我們得到了10,000個隨機變量的結果。如果您運行這段代碼，您可以看到null分布的形成是由於所觀察到的值堆疊在一起。
data(father.son,package="UsingR")
x <- father.son$fheight
n <- 100
population <- read.csv('https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/femaleControlsPopulation.csv')
library(rafalib)
nullplot(-5,5,1,30, xlab="Observed differences (grams)", ylab="Frequency")
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))
##if(i < 15) Sys.sleep(1) ##You can add this line to see values appear slowly
}
上圖相當於直方圖。從我們之前計算的null vector的直方圖中，我們可以看到obsdiff這樣大的值是比較少見的。這里需要記住的重要一點是，當我們通過計算案例來定義Pr(a)時，我們會學到，在某些情況下，數學給了我們Pr(a)的公式，省去了我們計算它們的麻煩。這種強大方法的一個例子使用了正態分布近似。
8、正態分布
我們上面看到的概率分布近似於自然界中很常見的一種:鍾形曲線，也稱為正態分布或高斯分布。當一組數字的直方圖近似於正態分布時，我們可以使用一個方便的數學公式來近似任意給定區間內的值或結果的比例:
雖然公式看起來嚇人,別擔心,你將永遠不會記憶,因為它是存儲在一個更方便的形式(如pnorm R,設置a為−∞,並以b作為參數)。這里µ和σ被稱為人口的平均值和標准偏差。如果這個正態逼近適用於我們的列表，那么列表的總體均值和方差可以用在上面的公式中。
一個例子是，當我們在上面注意到null 分布中上只有1.5%的null分布值在obsdiff之上。
總結，計算老鼠飲食差異的p值很簡單，對吧?但為什么我們還沒有完成呢?為了進行計算，我們從Jackson實驗室購買了所有可用的小鼠，並重復執行我們的實驗來定義零分布。然而，這不是我們在實踐中可以做到的。統計推理是一種數學理論，它允許你僅用樣本中的數據(即最初的24只老鼠)來近似這個結果。我們將在下面的部分中重點討論這個問題。
9、Setting the random seed
在繼續之前，我們簡要解釋一下以下重要的代碼行：
set.seed(1)
在本書中，我們使用隨機數生成器。這意味着，所給出的許多結果實際上可以隨機更改，包括問題的正確答案。確保結果不變的一種方法是設置R的隨機數生成種子。
練習random seed
dir <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/"
filename <- "femaleControlsPopulation.csv"
url <- paste0(dir, filename)
x <- unlist(read.csv(url))
這里x代表整個總體的權重。
10、Populations, Samples and Estimates
現在我們已經介紹了隨機變量、零分布和p值的概念，現在我們准備描述允許我們在實踐中計算p值的數學理論。我們還將學習置信區間和功率計算。
總體參數：
統計推斷的第一步是了解你感興趣的是什么。在小鼠體重的例子中，我們有兩個種群:控制飲食的雌鼠和高脂肪飲食的雌鼠，體重是感興趣的結果。我們認為這個總體是固定的，隨機性來自於抽樣。我們使用這個數據集作為示例的一個原因是，我們碰巧擁有這種類型的所有鼠標的權重。
dir <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/"
filename <- "mice_pheno.csv"
url <- paste0(dir, filename)
dat <- read.csv(url)
library(dplyr)
controlPopulation <- filter(dat,Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
hfPopulation <- filter(dat,Sex == "F" & Diet == "hf") %>% select(Bodyweight) %>% unlist
library(rafalib)
mypar(1,2)
hist(hfPopulation)
hist(controlPopulation)
我們可以用q-plot圖來證明分布是相對接近於非均勻分布的。后面將詳細討論，重要的是要知道它比較(y軸上的)數據和(x軸上的)理論分布。如果點落在恆等線上，則數據接近理論分布。
mypar(1,2)
qqnorm(hfPopulation)
qqline(hfPopulation)
qqnorm(controlPopulation)
qqline(controlPopulation)
兩個種群所有權重的分位數-分位數圖。
樣本量越大，對這種近似的缺點就越寬容。在下一節中，我們將看到對於這個特定的數據集，t分布甚至對於小到3的樣本容量也能很好地工作。
如果一組數的分布近似正態分布，那么隨機抽取一個數，它就會服從正態分布。If the distribution of a list of numbers is approximately normal, then if we pick a number at random from this distribution, it will follow a normal distribution.
dir <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/"
filename <- "mice_pheno.csv"
url <- paste0(dir, filename)
##The data has missing values.
##We remove them with the na.omit function
dat <- na.omit( read.csv(url))
y <- filter(dat, Sex=="M" & Diet=="chow") %>% select(Bodyweight) %>% unlist
avgs <- replicate(10000, mean( sample(y, 25)))
mypar(1,2)
hist(avgs)
qqnorm(avgs)
qqline(avgs)
11、實踐中的中心極限定理