Data Anaylsis for Life Science_I

2019-02-22 本文已影响2人城管大队哈队长

Introduction

这一个章节介绍了一些跟p-value以及置信区间相关的统计概念。

书中引用了某篇paper中的几段话：

“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.”

paper中，为了支持他们的结论，他们提供了下面的分析结果

“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.1 g) than in the normal dietfed mice (+ 0.2 ± 0.1 g; P < 0.001).”

这里的P-vaule以及 $\pm$ 指的是什么呢？在这一章节的后面，会慢慢介绍这些概念。不过首先，我们得先理解random variable（随机变量）这个概念。我们会使用一些数据来更好地理解这个概念。

原书中的老鼠数据下载链接已经失效了，新的数据在这里

dat <- read.csv("femaleMiceWeights.csv")

Our first look at data

这里的数据是用了24只老鼠，平均分成2组，一组用chow，一组用high fat (hf)饮食。在几周之后，研究人员对每只老鼠都统计了重量。

> 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

我们想知道的是，hf饮食是否会导致老鼠变得更重。因此我们可以对两组老鼠取平均值，然后看control和treatment组的均值差异。

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

obsdiff <- mean(treatment) - mean(control)
print(obsdiff)
## [1] 3.020833

我们可以发现hf处理的老鼠大约会重10%。但实际上，这并不严谨，因为这些均值实际上是随机变量。如果我们重复这次实验，还是各取12只做对照和处理，我们就会得到不同的均值差异。我们每次重复，我们就会得到不同的均值，不同的均值差异，我们把这类数据叫做随机变量。

Random Variables

在一部分，我们会继续探究随机变量这个概念。假设我们已经拥有了全部对照组的老鼠的数据，我们常称之为the population（总体）。通常来说，我们并不会得到总体的真实数据，但这里只是用来模拟而已。

population <- read.csv("femaleControlsPopulation.csv")
population <- unlist(population)

我们进行了三次随机抽样，来看看均值的变化。

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

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

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

可以发现均值是在不断变化的。我们可以不断地重复这一步，然后就会发现一些关于这个随机变量分布的有趣的事情。

The Null Hypothesis

现在让我们重新回到我们计算得到的均值差异 obsdiff 。我们会猜想，我们如何知道这个均值差异是由于饮食造成的呢？如果我们给24只老鼠都做同样地处理，是不是均值差异也那么大呢，还是会更小呢？统计学上把这个情况称之为Null hypothesis（原假设，虚无假设）。

原书中关于null hypothesis的定义很难描述。

The name “null” is used to remind us that we are acting as skeptics: we give credence to the possibility that there is no difference

因为我们已经有了总体的数据，我们可以多次地随机抽取24只老鼠，然后分成两组并给予它们同样的处理，看看这些多次抽样得到的均值。这里我们做了1000次循环。

> 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)
[1] 0.0127

在null向量中的值，我们称之为null distribution。后面会给出正式的定义。

我们可以发现在1000次模拟中，比 obsdiff 大的值只有很少的几次。所以我们得出了结论，即在没有处理差异的情况下，我们能观察到比 obsdiff 的值还要大的次数只占到1.27%。这就是我们所谓的P-vaule，后面会给出更正式的定义。

Distributions

分布可以认为是一种大量数据的简单描述。例如，你测量了人类所有男性的高度，那么该如何展现这个数据呢？

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

最简单的方法就是简单地抽取一些数据

round(sample(x,10),1)

Cumulative Distribution Function（累积分布函数）

简单地抽取一些数据可能并不完整，所以我们就可以引入累积分布函数来展现数据的分布，即在我们计算的所有数据中，数值小于a的比例。
$F(a)=Pr(x\le a)$
上述公式就是所谓的累积分布函数（CDF）。当CDF来自数据，则称之为empirical CDF （ECDF）。我们可以画出累积分布图

hegihtecdf <- ecdf(x)
plot(hegihtecdf)

原书中的代码我不太能理解，所以没有照搬原书代码，但图画出来是一样的。

ecdf这个函数比较特殊，其返回的是一个函数，并不是像平常R函数一样返回数值。

Histograms

这里画的应该是频数分布图

hist(x)

Probability Distribution

我们也可以描述随机变量的概率分布情况。
$Pr(a \le X \le b)=F(b)-F(a)$
注意这里的X是大小，代表的是随机变量的概率分布，区分于小写的x，即随机变量。

举个例子，我们可以描绘出之前提到的，基于零假设成立的情况下，两组均值差异的分布图。

n <- 100
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
}

也可以换个图来看，可以看到跟先前我们计算出来的obsdiff一样大的值很少。

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

对于离散型的随机变量而言，展现的方式有：累积分布函数和概率分布函数

对于连续型的随机变量而言，展现的方式有：累积分布函数和概率密度函数

Normal Distribution

我们之前所见的概率分布函数是一种在自然界非常常见的分布：钟形曲线，也叫正态分布或高斯分布。当数值的分布近似于正态分布之后，我们就可以用积分公式得到在任意给定区间内的概率
$\mbox{Pr}(a < x < b) = \int_a^b \frac{1}{\sqrt{2\pi\sigma^2}} \exp{\left( \frac{-(x-\mu)^2}{2 \sigma^2} \right)} \, dx$
之前的例子是做了1000次循环，然后计算了比obsdiff大的值所占的比例。这里我们认为我们列表中的数值逼近于正态分布，就可以用函数计算了

1 - pnorm(obsdiff,mean(null),sd(null))
## [1] 0.01391929

可以发现结果是差不多的。

Summary

我们之前的多次循环等同于进行了一遍又一遍地实验，这是很耗时耗力的。但统计推断（Statistical Inference ）可以只根据样本信息，比如最开始的24只老鼠，让你近似做到这一点。

Setting the random seed

因为我们会用到多次的抽样操作，所以有必要设置种子，来让你的结果不会改变。

set.seed(1)

Exercises

x <- unlist(read.csv("femaleControlsPopulation.csv"))

设置种子数，从总体中每次抽取5个样本并计算样本均值，重复1000次。在这1000次中，大约有多少样本均值是偏离总体均值1 ounce的。如果重复了10000次呢。
```
n <- 1000
null <- vector("numeric",n)
for (i in 1:1000){
  null[[i]] <- mean(sample(x,5))
}

mean(abs(null - mean(x)) > 1)
```
大约是50%是偏离总体均值1 ounce的，重复1000次和10000次是差不多的。
如果设置样本数量为50呢，偏离总体1 ounce大约有多少呢。

就只有大约1.8%了。
画出条形图来查看，发现样本数量为50的，图形更加尖。
```
hist(null)
```
对于样本数量为50的均值分布，均值为23到25之间，大约有多少呢。同样地，对于均值为23.9，方差为0.43的正态分布，大约有多少呢。
```
mean(null < 25) - mean(null < 23)

pnorm(25,23.9,0.43) - pnorm(23,23.9,0.43)
```

Populations, Samples and Estimates

Population parameters

我们恰好拥有总体，可以计算一些总体的参数

dat <- read.cav("mice_pheno.csv)

library(dplyr)
controlPopulation <- filter(dat,Sex == "F" & Diet == "chow") %>%
select(Bodyweight) %>% unlist
length(controlPopulation)
## [1] 225

得到的是对照组的总体，把上述得到的225个数值分别记为 $x_1,\dots,x_m$

hfPopulation <- filter(dat,Sex == "F" & Diet == "hf") %>%
select(Bodyweight) %>% unlist
length(hfPopulation)
## [1] 200

得到的是处理组的总体，把上述得到的200个分别记为 $y_1,\dots,y_n$ .

就可以计算总体的一些参数了，例如均值和方差(variance )
$\mu_X = \frac{1}{m}\sum_{i=1}^m x_i \mbox{ and } \mu_Y = \frac{1}{n} \sum_{i=1}^n y_i$

$\sigma_X^2 = \frac{1}{m}\sum_{i=1}^m (x_i-\mu_X)^2 \mbox{ and } \sigma_Y^2 = \frac{1}{n} \sum_{i=1}^n (y_i-\mu_Y)^2$

然后我们就可以提出问题了： $\mu_Y - \mu_X = 0$ ？

Sample estimates

先前我们对于每个总体都获得了样本量为12的样本。我们把这些数据写作大写，说明这些数据是随机的。所以样本是 $X_1,\dots,X_M$ 和 $Y_1,\dots,Y_N$ 。 $N=M=12$ 。而上面的总体数值，其不是随机产生的，所以是小写。

因为我们想要知道的是 $\mu_Y - \mu_X$ 是否等于0，所以我们考虑样本均值差是否为0。
$\bar{X}=\frac{1}{M} \sum_{i=1}^M X_i \mbox{ and }\bar{Y}=\frac{1}{N} \sum_{i=1}^N Y_i.$
注意均值之差实际上也是个随机变量。先前，我们已经了解过随机变量了。通过不断地从原始分布中抽样。但实际上，在实践过程中，这是不可能存在的，那相当于你不断地重复买24只老鼠。在这里，引入了一个数学理论。在数学理论上，将 $\bar{X}$ 和 $\mu_X$ 以及 $\bar{Y}$ 和 $\mu_{X}$ 关联起来，反过来帮助我们理解 $\bar{Y}-\bar{X}$ and $\mu_Y - \mu_X$ 。接下来也会理解中心极限定理 (the Central Limit Theorem) 和t分布。

Exercises

library(tidyverse)
library(rafalib)

dat <- read.csv("mice_pheno.csv")
dat <- na.omit( dat )

# 1
x <- filter(dat, Sex == "M" & Diet == "chow") %>%
  select(Bodyweight) %>% unlist()
mean(x)

# 2
popsd(x)

# 3
set.seed(1)
value_x <- mean(sample(x,25))

# 4
y <- filter(dat, Sex == "M" & Diet == "hf") %>%
  select(Bodyweight) %>% unlist()
mean(y)

# 5
popsd(y)

# 6
set.seed(1)
value_y <- mean(sample(y,25))

# 7
abs(value_y - value_x)
abs(mean(y)-mean(x))

# 9 
# 相比于M，Female 样本估计得到的均值差异更接近于总体的均值差异。可能的原因是
# C：female的样本更大，更能模拟总体。

Data Anaylsis for Life Science_I

Introduction

Our first look at data

Random Variables

The Null Hypothesis

Distributions

Cumulative Distribution Function（累积分布函数）

Histograms

Probability Distribution

Normal Distribution

Summary

Setting the random seed

Exercises

Populations, Samples and Estimates

Population parameters

Sample estimates

Exercises

猜你喜欢

热点阅读