假设检验 以及 qPCR数据处理应用

2020-03-28  本文已影响0人  小折线

前言1：
首先 向肖恩的好读者，好伙伴们道一声感谢，感谢一路的陪伴，最近有数位读者来询问问题，让我很开心。在此，也需要向大家道歉，大家在微信公众号的留言，只有48h内可以回复，而我是个学生，不能保证固定时间查看公众号。而且微信公众号没有提醒，有所怠慢请见谅咯。

前言2：
2020是新的一年，然而，这新年伊始充满了坎坷，充满了让大家措手不及的意外，但我相信大家对生活的热情不减，生活还要继续，既然如此，何不重整行囊，砥砺前行。（毕竟我已在家葛优多日，是时候学习啦）

那么进入正题吧！

今天我们讲讲假设检验这个东西，并且会以qPCR数据处理为例子，应用假设检验，获得传说中小于0.05就很了不起就P值。

文末还有大大的彩蛋哦！！！

假设检验

假设检验这个东西，我们顾名思义，就是检验某一个假设，所以，这就很自然地引出两个重点，欸对，两个重点就是：

  1. 检验 （检验的方法）
  2. 假设 （检验的目的）

我们可以这么理解student t检验这个东西：

我有两组数据，A和B。我想知道A组数据和B组数据的均值是不是有显著性差异。 那么我们要先有一个假设，然后检验这个假设成不成立

一般来说，我们的假设是，A，B两组没有差异，然后计算在A，B没有差异的情况下，出现我们观察到的数据的概率。如果该概率<0.05,那么我们就说，这事发生的概率太低了，这个假设不靠谱，我们就否定原有的假设，进而相信A组数据和B组数据在统计上有显著性差异。

这里有两个统计学的术语：

除了T检验，还有别的检验，不同的检验都有对应的假设和检验方法，比如F-检验，卡方检验，秩和检验等，适用于不同的情况，有不同的目的，使用前要搞清楚。

Student T Test（T检验）

最常用的就是T检验，我们已经知道了他的整体逻辑了，但是我们还不清楚，他如何计算的，理解了计算过程，可以让我们更好地理解和解读p值，不去错误地使用p值，不迷信p值。

在理解计算原理之前，我们举个栗子，思考一下别人的人生：

B站的罗翔有一个CP，叫张三。

张三是个果农, 种了两种苹果树，一种叫A，一种叫B，他想知道哪种苹果树结的果子大。他在果子成熟时，分别摘了A苹果100个，B苹果100个，发现A苹果平均300g，B平均305g。他于是下结论：B苹果树结的苹果比A大。于是计划着减少A的种植数量。这时候，他的脑海中出现了一个声音 “老张你这不太严谨哦”。老张吓了一跳，"见鬼啦，谁在说话",一个富有磁性的嗓音说：“我是来帮你的，如果想要苹果长的大，还请老张思考以下问题:”

  • A苹果300g，B苹果305，你觉得B比A大，那么如果A平均300g,B平均301g呢，如果B苹果300.5g呢，你还觉得B比A大呢?

  • 会不会出现以下情况，你两次都从同一棵A苹果树摘苹果，有没有可能，第一次300g，第二次305。就算同一棵树你两次摘苹果，结果都可能不一样的，两次摘苹果，平均值差0.5g，差5g，差10g都是完全可能发生的，那平均值差多大才算大呢

  • 最后，请思考以下两个场景，

    1. 你从AB两种树，分别摘10个苹果，B重量平均比A大5g
    2. 你从AB两种树，分别摘100个苹果，B重量平均比A大5g

    你觉得哪一种场景下，B苹果树的苹果重于A的概率更大。

老张闭上眼睛，认真思考着以上问题（读者可以跟着他一起思考，助他一臂之力）

没想到老张是个被苹果耽误了的数学家。他花了一周画了两张图。

# 不感兴趣的同学不需要看
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
# fig 1
norm1 = stats.norm(300,3)
norm2 = stats.norm(310,3)
x = np.linspace(250,350,100)
y1 = norm1.pdf(x)
y2 = norm2.pdf(x)
fig, axs = plt.subplots(1,2)
ax1 = axs[0]
ax1.plot(x,y1, color='b')
ax1.plot(x,y2, color='r')
ax1.set_ylim((0,0.15))
ax1.axvline(x=300, color='b', alpha=0.5)
ax1.axvline(x=310, color='r', alpha=0.5)
# fig 2
norm1 = stats.norm(300,30)
norm2 = stats.norm(310,30)
y1 = norm1.pdf(x)
y2 = norm2.pdf(x)
ax2 = axs[1]
ax2.set_ylim((0,0.15))
ax2.plot(x,y1, color='b')
ax2.plot(x,y2, color='r')
ax2.axvline(x=300, color='b', alpha=0.5)
ax2.axvline(x=310, color='r', alpha=0.5)
fig1.png

老张点了根烟，叹息道，这一切都是抽样造成的，如果我能收集世上所有A苹果和B苹果，就可以求出准确的均值，我就知道哪个苹果大了。但是我做不到，做不到就只能抽样。抽样只能近似真实的均值，但总存在误差。

比如上面的两张图，蓝色是抽样的A苹果不同重量的频数，红色是B苹果。左右两张图，同样是均值差10g（均值以竖直直线表示），但是左图中，红色曲线下面积与蓝色曲线下面积很少重叠，也就是说，红色整体上确实大于蓝色。而右边的图，虽然红色的均值大于蓝色10g，但是这两条曲线基本重合，这次抽样B比A大10g，很有可能是误差。下次再试一次，可能A就比B大了。

老张猛吸一口烟，转念又想，即使是右图的情况，如果我是统计了咱们镇所有的AB苹果，得出这样的结果，是不是还是可以说B比A平均大10g呢，这总比我统计100个苹果来的准吧！

老张熄灭了手头的烟，缓缓的说：

“失去苹果，失去很多；失去统计学，失去一切。”

面对浩瀚无垠的数学宇宙，张三开始了他的统计之旅：

他提出了以下思路：

  1. 我们的目的不是比较被抽样的苹果均值大小，而是想利用抽样来比较两种苹果树真实的苹果均值的大小
  2. 我们对苹果的真实均值的大小只能估计，我们的估计不应该是一个值，而是一个范围，或者一个分布，一个均值的概率分布
  3. 我们对真实均值的概率分布估计必须考虑数据的离散程度（即方差）和抽样次数。方差越小（苹果重量波动小），抽样次数多，那么估计就越精准

于是乎，老张不知道怎么回事（我也不知道怎么回事，请知道这个公式怎么来的同学私戳我），推导出了第一个公式

se(\overline{x} ) = \sqrt{\frac{\sigma^2}{n}}\approx \sqrt{\frac{s^2}{n}}

x表示苹果的重量，\overline{x} 表示苹果的均值，
\mu 表示苹果的真实均值， n为样本数
\sigma为x的真实标准差，s为x的样本标准差

se(\overline{x} ) 表示\overline{x} 与\mu 之间的误差的标准差，学名叫标准误，standard error （se）

该公式是说，如果我对一颗苹果树抽样，取n个苹果，那么我们可以认为苹果的样本均值概率应该服从
以\mu为中心,标准差为\sqrt{\frac{s^2}{n}}的正态分布

注意： 这里描述的是采样均值的分布，不是采样个体的分布，即你采样很多次，每次计算出的均值放在一起，是这样一个正态分布

这里可能难以理解，于是老张又画图说明

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
norm1 = stats.norm(300, 10)
s = norm1.rvs(30)
x = np.linspace(270,330,100)
mean = s.mean()
std = s.std()
mean_std = np.sqrt(std ** 2 / 30)
y2 = stats.norm(mean, mean_std).pdf(x)
ax = plt.subplot()
ax.hist(s,bins=10, density=True)
ax.plot(x,y2, color='r')
fig2.png

这张图这么理解，这里很重要，
我是上帝，我知道苹果树的苹果真实均重300g，标准差10g，我从一种树上摘了30个苹果，对苹果的取样应该，其均值应该是以300为中心，标准差为\sqrt{\frac{\sigma^2}{n}} = 1.8的正态分布。

张三不是上帝，张三只能根据采样，他发现样本均值为298，标准差为7.78，（图中蓝色部分）。但是他对真实的均值和标准差一无所知，所以聪明的张三，反其道而行，他说，既然能用真实的均值和方差推测样本的均值，是不是可以用样本的均值和方差推测真实均值，套用上面的公式，推测苹果的真实均重概率应该是图中红线那样，以298为中心，标准差为\sqrt{\frac{s^2}{n}} = 1.42， 这就是对真实均值力所能及最好的判断

这时候，张三转念又想，怎么比较A和B的差异呢，眼珠子那么一转，有了！

我们既然比较AB之间的差距，为什么不去计算 A均值-B均值呢，
张三说要有X，我们就有了X
X = \overline{A}-\overline{B}

这是，我们只要检验 X = 0 是否成立就行了，我们的零假设是AB无差异，即X=0， 如果零假设成立，那么X的分布应该是一个以0为中心的分布，该分布的方差，应该是 \overline{A}方差 +\overline{B}方差 ，即：

se( \overline{A} \pm \overline{B}) = \sqrt{ \frac {s_A^2}{n_A} + \frac{s_B^2}{n_B}}

有了这样的分布，我们就可以描述每次抽样X在某个范围的概率

张三画了第三张图

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
norm1 = stats.norm()
x = np.linspace(-4,4,100)
y1 = norm1.pdf(x)
ax = plt.subplot()
ax.plot(x,y1, color='b')
l1 = norm1.ppf(0.025)
l2 = norm1.ppf(1-0.025)
ax.fill_between(x,0,y1,
                (x < l1) | (x> l2), color='darkred')
fig3.png

该曲线下面积占比其实就是零假设下的事件发生的概率。

P值计算的是零假设下的发生某事件或者更为极端事件的概率

如我取样后计算得到X=2， 按照正态分布可以计算得，从2到正无穷，曲线下面积占比为约0.025（右边红色部分），那么A比B大2或者大更多的概率是2.5%，这就是P值。2.5%这个概率够小，我们觉得零假设应该是在扯淡，于是拒绝零假设，选择备选假设，即认为AB有显著差异

上面其实是算单尾的情况，即比较A比B大，或者A比B小的情况，会先假定一个方向。而日常常用的是双尾t-检验，他不假设AB哪个大哪个小，所以如果算出X=2， 他会计算两边，(-\infty,-2),(2,\infty)的曲线下面积，故数值上P值=0.05，是单尾的2倍

一般我们以5%作为阈值，只有在P<0.05时，我们才认为，零假设下，出现这么大的X的概率太小了，才会舍弃零假设，选择AB确实有差异这个备选假设。
因为双尾不假设AB哪个大哪个小，p值是单尾的两倍，所以更为严格。我们一般都选双尾

重要说明：前面其实我骗了你们，刚才我讲的，严格来说，不叫t-test，t-test所用的不是正态分布，而是t-分布，t-分布类似于正态分布，比正态分布要胖一些，而且受自由度控制（采样量-1）。一般来说采样量 > 30时，用正态分布，采样量 < 30时，用 t - 分布。

qPCR数据应用

说了这么多，我们来应用一下，比如我有基因A，我一顿操作处理了细胞，想通过qPCR看看A的表达变了没。
比如 我来编个数据，

条件 重复1 重复2 重复3
未处理 1.5 1.5 1.6
处理 1.6 1.7 2

那你说这个A的表达是变了还是没变？

我也不知道，但我们可以假设A没变，看看我们获得这样的数据的概率，如果概率<0.05, 我们就舍弃这个假设。

在Excel中，可以直接算t-test,如图：


t-test.png

尽管上下两次计算的时候，他们的均值差大小相近，但由于第二次样本多，所以我们对真实均值的估计分布会精确一些，即正态分布更加修长苗条，他们同样的均值差情况下，第二种情况就落在正态分布更外侧的地方，曲线下面积更小。

福利时间

很感谢，大家看到这里，还记得之前我出过一个qPCR处理的软件吗，这次推出了进化版，可以计算p值，程序会自动计算每个基因处理组与对照组的双尾t检验，还是一键全自动，熟悉的味道，更强的功能。同名公号后台回复qPCR领取。

zzzz.gif

2020-03-28

上一篇下一篇

猜你喜欢

热点阅读