不止是线性回归

序言

很高兴收到大家的互动，虽然是系统自动回复，但正好今天打开公众号，想写的是，这篇三个月多就一直在纠结、探索和求证的事情。

今天来想讲一下差异分析。说实话，我虽然从事了10年生物信息，博士期间主要课题是构建哺乳动物转录组数据库，严格意义来讲，其中的分析就是跑流程，为了给用户提供广泛的数据资源以供用户自由筛选。博士后期的多维数据整合分析，倒是学会了各维组学数据的分析，在结果的归纳处理，统计学方面的学习和运用还是有所欠缺。

第一个课题做得页面很炫酷，呕心沥血地用了很多高大上的技术，导致，我时常在想，什么是生物信息，为什么结合实验的文章看上去毫无技术含量，为什么结合了实验就能发表好文章？至少博士期间我是这么愚昧的认为的。我还时常跟师姐讨论自己内心中对差异表达分析的不屑。因为在之前的我眼里，例如基因表达的差异分析，就是当时的主流软件主流流程（tophat,cufflink)一甩，几个实验条件一设计好，就出结果了。生物信息费老大劲，纯信息大于5的杂志算好杂志。当时就觉得，差异分析，毫无技术含量。然而，现在真的从事差异分析，才发现并非如此。

差异（表达/甲基化...)分析/也或者各类qtl分析...

做这类分析的方法很多，之前说过，广泛使用的有rank sum检验、线性回归及其各种衍生算法、t检验等等，尽管可以使用的软件和R包层出不穷，都是使用的工具，故万变不离其宗。在我做这类分析的时候，我就问了自己很多问题，rank sum和t test可以出很多结果，百来个到千把个结果。linear regression出的结果是需要技术，F检验出2个结果。那么问题来了，什么结果才是可信的呢？纵观各类杂志，linear regression找差异普遍发的杂志很好，当然这并不是终究原因，最终原因当然是结果好！可信！根据对算法的理解，可以说rank sum和t test都是很less power的，找出几千几百个，其中90%都是假阳性的结果，在此基础上产生什么揭秘疾病的分子机制类的结论，是不是有些忽悠，作为作者，就问自己的内心，你信你的结果吗？如果你自己可以在内心中确信，那就是好结果。（后续详述）

我想学分析，我想学R，R好学吗？

一言以蔽之，不用学R，R只是工具，拿着就用。就像线性回归，你知道统计学算法原理，用自变量回归因变量，拟合出一条直线，有了斜率，再做F, t或者r检验，看看线性关系是否显著，再看看模型是否拟合得好，考虑的因子是否充分，残差是否服从正态分布。嗯，这是一般思维。那么在R里面，就是一句话lm(y~x)。需要学R吗？那就，我现在需要用线性回归，我直接百度或谷歌一下，线性回归R,就有成千上万个中英文帖子告诉你，lm一下，结果立现。R就是个黑盒子，学一下变量、数组、数据框之类的数据结构就可以了，其他的需要什么功能，实时百度。难道非得把所有函数名字背下来，才够格对外宣称，我会R？就像使用遥控器一样，难道还需要看使用说明？

说说多元线性回归

举个例子，例如年龄和某个疾病有关，我设置x为年龄，y为疾病的生理指标（假设这里是服从正态分布的），我lm(y~x)，那所得的斜率就代表他们相关？答案当然是否定的。当我们加入另一个主导因子c的时候，必然会发现c的出现，削弱了年龄的斜率，甚至能让年龄和生理指标的线性关系不再显著。这就是为什么，各类差异分析或各类qtl找相关类分析，都需要去掉confounder/covariate了。因为必须去掉相关因子的影响，最后才是真正你关心的例如年龄对某生理指标的影响。

什么样的结果才是好结果？相关类分析中，我们应该关心的是什么？

不用想的回复就是，每次都显著，无论怎么折腾。

网上帖子都会讲线性回归。然后附加很多其他需要关注的东西。例如残差。如果模型拟合得好，残差就应该服从正态分布（因为就相当于是随机误差），说明模型拟合完全。然而，例如我只想找某个探针的甲基化水平和某生理指标是否相关，我需要把所有能解释因变量的因子都抠个遍，直到残差只剩下随机误差为止吗？那么我的模型解释的比率R2，岂不是能高达80%。这才不是我们关心的结论。我们关心的是，是否相关！！！

肯定是当我把所有会影响我所关心的因子和因变量的关系的所有其他因子抠出来，直到其不再削弱关心的因子和因变量直接的斜率显著性关系为止。

那么，我们要怎么去知道，是否抠干净了呢？可以尝试作一个很重要且被人们熟知的qq图。众所周知，qq图是用来简单看看实际分布和理论分布是否一致的。设想，如果是随机分布的，最终得到的概率p值，应该是均匀地分布在[0,1]区间，那么，如果所有的影响两个元素相关性的中间干扰变量我都已经排除完了，最后的p值真正能证明我的结果，那么绝大多数其他不显著的p就应该是服从均匀分布的。这个图可以非常轻松的用R包qqman来实现。

倘若图形是这样，最明显的上面这个点的p值就远远显著过其他点，而这个点的结果，已经在其他人的两个实验中都有重复发。其他点都是服从均匀分布，落在这个直线上，可以量化这个信息，计算膨胀系数，就是对两个分布的统计值的比较，如果膨胀系数约接近于1，就说明，理论分布越符合真实分布。（这是判断结果是否可信的基准）

那么更一般的结果是什么？很多能找几百几千个的结果，拿出p值画qqplot大约应该是这样的

严重膨胀，这就说明confounder还是没有找齐全，甲基化和性状之间的p值还是会随着confounder的参入而改变。结果不可信。

真正好的结果，是各类统计图都能从方方面面去证明的。直到文章原作者本人都深信。甚至构建模型好几个月，显著的永远都是显著的。我想这就逼近真实了。

就一个线性回归，又不高大上，有什么了不起？

曾经旁听朋友跟一个搞算法的朋友聊天。虽然个人一听“搞算法”三个字，就觉得这人牛逼哄哄，至少我搞不起，因为我觉得搞算法就得像真正的算法大师一样，修正算法中的不足。而不是，在这套数据中，我用了遗传算法，那套数据中，我用了广义逆高斯模型，贝叶斯网络，svm，听上去，我就牛逼哄哄。实际上，结果不可信，算法reviewer看不懂，做对做错了还是看不懂，对科研是否有帮助，自己也不知道，只剩下，这些炫酷的名字，看上去让自己牛逼哄哄。

周志华先生的一本很好的西瓜书，把机器学习讲得浅显易懂，可谓大师。西瓜书的开场是介绍线性回归。因为线性回归其实是svm等等各种算法的基础。所以属线性回归运用最广、最基础。能解决的问题最多。倘若同样能解决问题，尚且每个名字牛逼的算法背后都有一定局限，为何非要不屑线性回归。

对模型的改进方法

差异分析的难点就在这里了，怎么才能构建出好的模型。当然，我觉得最有帮助的就是找相关这方面的医生聊天，医生能一眼看出什么样的生理指标和你要研究的性状相关。但是，记住，真相永远是可以相互、反复印证的！

医生的先验知识源自于一代代的长时间观察、统计、归纳

倘若你用统计学方法也能找出来，这时候你可以和医生去对特征，你会发现，跟你找出来的几乎相似。合情合理，和看上去通篇都是理所当然的结论，就足够好。

在筛选模型特征的时候，可以仔细考量每个自变量，是随机效应，还是固定效应。可以按具体情况适当修改线性模型，例如考虑使用线性混合模型。

可以用向后或者向前，个人喜欢向后渐少特征的方法，用方差分析对增加或减少特征后的模型的变化的显著性进行比较下结论。

比较快捷的方法，推荐斯坦福大师团队推出的LASSO方法进行特征筛选。

再回归线性回归的基础思考，画出因变量的概率密度图，线性回归要求因变量符合正态分布，这往往也是造成膨胀的因素。所以通常在用基因表达量作为因变量的时候，往往会做rank normal transfer以构建服从正态分布的统计量。

关于sva方法，其实推断sva的时候，矩阵的构建也是需要仁者见仁、智者见智的，但是，三个月实践下来，觉得sva并不是很好的有助于模型提升。另外，sva的引入，会导致审稿的时候难以解释清楚。

另外需要注意的是，在做各种转换的时候，会损失信息量，会让模型和结果的说服力降低，听者会觉得改变了所要找的本质，会less power。也就是能不变换，就避免变换。

由此可见，同样是R中的一句lm，也或者同样是一个matrixEQTL，不同的努力，得出的结果和可信程度以及发表的文章质量是完全不一样的。这就是差异数据分析的难点。