limma 算法总结

我这个人完全没办法在不了解算法的情况下用一个包，那样我会疯，可是网上居然完全找不到limma算法的讲解资料，全是介绍怎么用的。

可能我太初级了，这么简单的算法大概不值得一提……
于是搞了一天原文，和东拼西凑的统计学知识，把现在了解的先总结一下：
limma是做差异表达的包，其算法核心在两个function上：
lmfit（）以及eBays（）

lmfit（）就是multiple linear regression

假设我们有基因表达矩阵$Y=[y_1,y_2,...y_{100}]$,其中包括50个样本和100个基因，样本包括30个病人和20个健康的
这样病人vs健康可以写成一个长度为50的0-1向量X=[11111111….0000000], 1表示病人，0表示healthy control

那么算法的大体思路是：

1. 假设数据质量良好，我们只想知道病人和健康人群的差异表达，对于每一个gene i，构建 一个线性模型
   $y_i=ax+b$，然后利用最小二乘，求出这条线的a和b；这就是lmfit的作用。最小二乘法就是最小化实际的$y_i$和拟合直线上的${\hat{y}}_i$的绝对距离。

2. 由于只有一个变量，那么这就是个简单的线性回归。（这个时候，机器学习和统计学的区别就出来了：如果是机器学习，找到一条线，就要开始交叉验证，在测试集上验证模型的鲁棒性，然后选出最好的模型来。但这些都是基于我们有大量的样本允许我们这么折腾。而很多时候我们就没几个样本，所以生信啊，大多还是基于统计模型。另外，'交叉验证，测试集’什么的主要是为了说明模型的鲁棒性，目的是为了拿这个稳健的做预测，你说它这有什么理论依据，只能拿测试集结果的‘准确率’来说明模型有效。） 扯远了，由于这是个统计模型，模型的有效性需要各种假设检验来证明。我们的目的是找差异表达基因，拟合了线性模型，怎么看我们的x和yi有多相关呢？就需要用到$R^2$
   $R^2=(var(y_i)-var(y_{fit}))/var(y_i)$,
   yi的方差就是说，不考虑任何x的因素，单看yi上的样本离均值有多远，yfit就是考虑了x之后，我们的样本离拟合线有多远， 如果x和yi很相关，那么加上x的因素之后，yi的方差就会变小，$R^2$就会很大 。（p.s.这里机器学习和统计的不同在于，机器学习才不在乎什么相关不相关，它一般只在乎拟合的好不好，能不能用这个线去预测别的人有没有病）

3. 然而，统计模型的麻烦之处就在于，你不仅要找到有多相关（$R^2$），还要验证你找到的相关性有多显著。这里用F分布来求出p-value来说明,模型的显著性。（为什么是F分布？）

4. 注意F分布的p-value是为了测试模型的显著性，我们这里只有一个变量，所以模型的显著性也就是x的显著性了。但是如果有两个或以上样本，F test是为了测试单纯的yi的方差，和用所有样本x1,x2…拟合之后的$y_{fit}$差异的显著性。

   t test: This test checks the significance of individual regression coefficients.
   F test: This test can be used to simultaneously check the significance of a number of regression coefficients.

5. 可是，事实上我们可能并不关心模型的显著性，我们只是想知道哪些基因和某一个疾病（特征/自变量）相关。尤其是，对于limma，如果我们的数据存在‘batch effect’，我们还需要把batch也作为自变量输入进去，对于batch这一列自变量，我们根本就不关心它拟合的好不好，甚至还有点不想让它存在。因此我们真正需要的是对我们感兴趣的单个自变量进行的t-test检验，

6. 对某一个自变量x做t-test，就是测试，如果我们拿掉x（线性回归参数a=0），和不拿掉x（原来的拟合结果）差异的显著性，这些结果对应了imma里的t-value和p.value

7. 多重假设检验: 然而还没有结束，因为我们不止对一个gene做了回归，我们对全部gene都做了回归，gene i的t-test的pvalue表示了结果的假阳性率，可是所有的gene都计算，假阳性总和就会变很高，结果就不可靠了。因此我们必须要进行 multiple hypothesis test。以BH方法为例，比如100个gene，ttest之后我们觉得结果里面有假阳性，所以就对这100个gene按照pvalue从小到大排序，然后计算一个新的qvalue：
   $q_i=p{v}_i\ast (m/k)$ 其中m是总共检测数量（m=100），k是排名
   这里我用自己的数据画了一个pvalue和adjusted pvalue的关系（m=10000），可以看到这个凸函数，放大了pvalue中较小的值
   

eBays（）就是用经验贝叶斯调整了t-test中方差的部分

先贴原文：

对比一下经典的t-test：

可以发现这里，就是用后验${\tilde{s}}_j$去替代了原始的$s_j$.这才是limma的主要贡献，以上都是基本的多元线性回归。可是为什么要这样做？

没时间了，过两天补上。。

Reference

http://reliawiki.org/index.php/Multiple_Linear_Regression_Analysis#:~:text=Test%20on%20Individual%20Regression%20Coefficients%20(t%20Test),-The%20t%5C%2C%5C!&text=test%20is%20used%20to%20check,may%20make%20the%20model%20worse.

https://bioconductor.statistik.tu-dortmund.de/packages/2.11/bioc/vignettes/limma/inst/doc/usersguide.pdf