对于R语言中的model.matrix的一些简单理解

首先，需要简单的申明一下，model.matrix是R语言自身的包，不是bioconductor包里的部分。

所以在理解上，不需要带上生物信息学的知识，也无需生搬硬套bioconductor里的其他函数的功能。

然后，我对于这个函数的理解还是有不足的，特别是factor到底是怎么处理的。后来去看了源码，但是源码有点长，看不太懂。

正式开始吧，

这个函数，他是应用在线性回归时自动生成设计矩阵的，那么肯定大家都好奇到底什么是设计矩阵。

我搬用“https://genomicsclass./book/pages/expressing_design_formula.html”上有人写好的一篇文章里的公式，

如果我们要对形如这样的线性方程进行回归，也就是计算β0和β1的值，那么：

改写成这样后，这个X就是我们的设计矩阵。

所以例如x  = c(1,2,3,4,5,6)的时候 我们调用model.matrix(~x)会生成

这个函数的含义是以x为自变量的一阶方程的设计矩阵

形如这样的设计矩阵就比较好理解了，刚好是同我们的定义是吻合的。

实际上！我们在limma学习过程中，并不是真的在进行线性回归，我们在做的工作叫做差异分析！

TypeScript

我觉得虽然在公式上，某一些简单的条件下这两种事情的做法是一模一样的，但是还是要弄清楚。
因为我再学习的时候，看的资料一会儿说线性回归，一会儿说差异分析，弄得我根本不知道到底目的是什么。
而且狭义上的线性回归根本不是差异分析，本身在计算过程中的意义就不一样。--------来自我的吐槽，希望大家不要搞混，下面有解释

可能大家到这里还是不太理解，这个矩阵好像和我们在limma里面生成的完全不一样。

是的！的确这也说明了我们不是在搞什么线性回归，所以不是直接拿这一套矩阵就能用的。

首先我要先解释，在生物信息学的limma包中，为什么要转化为因子（factor）呢？

在解释之前我想先把差异分析说明白，差异分析就是将数据按照一定的要求进行分组，然后比较不同组之间的差异。

这个分组的方式称之为分组矩阵，实际上就是我们这里用到的设计矩阵。

从本质上说，我们这个分组矩阵叫做设计矩阵太勉强了，而且很容易把大家搞混。

而比较的方式用比较矩阵进行定义（比如说哪一些数据更重要一些，哪一些不重要的权重划分以及怎么比较等等功能）。

好弄明白差异分析之后，因为我们在用limma做差异分析的时候，需要分组，而分组的自变量往往不是数字，比如说sCLLex数据里，我们是对发病还是不发病作为分组的自变量。

此时，我们转为因子，就是将分组的自变量从一些字符类型转化为了多维坐标上的一个点。

例如factor(c('我们','你们','他们'))就是把“我们”作为(1,0,0)“ 你们”作为(0,1,0)“ 他们”作为(0,0,1)

而且我们需要的是分组，所以必须用一个统一的数字去表示，这里是用的1，那比如我们统一用2去表示当然也可以的，比如说

(2,0,0) (0,2,0) (0,0,2)在使用上是不会影响真正的结果本质的，当然可能数值上有些许差别。

（实际上并不是这样，但是这样能更好的理解，实际上是吧因子中每一个值都变成了一个自变量，比如说factor(c('我们','你们','他们'))其实是有三个自变量

这样的话设计矩阵就是这三个自变量进行回归时候的设计矩阵了，如果这个值出现了，那么这个自变量为1，最好回头去理解这个因子，其实这只是我自己的理解，可能不正确，只是现在凑巧这种理解方式能够解释为什么这个函数会生成这样的结果罢了）

例如我们的输入为x = factor(c(1,2,1,2))时，调用model.matrix(~x)的到的结果是

这种做法很好的进行了分组，而且消除了不同组之间（比如实际上我们的自变量是组1和组2而不是数字1，2）而带来的错误。

在理解上，我们回到公式，在组1中的数据进行回归时，带入公式得到y1=β0 因为组一的自变量等于说是0 因此得到的β0就是组1的平均值

组2的数据进行回归时，带入公式得到y2=β0+β1，所以β1=y2-β0，所以β1表示的是组2的y值减去组1的y值平均值后再平均的结果。

那么有趣的是刚好满足我们的求组间差异的要求，这里的β1就是我们的组间差异了。所以根本用不到比较矩阵，直接进行线性回归就行了。

但是这种方式的差异分析显然是很简单的，比如权重啊什么的通通都没法考虑了，完全是平均数相减。

所以一定要记住差异分析和线性回归不是一回事，只不过我们凑巧有相似的部分罢了。

接下来是说我们在limma使用时用的~0+factor()的方式，这种方式我到现在还是不理解为什么会得到按组分离的矩阵，我现在的理解就会这个函数对这种情况 做了特殊的处理。

也就是按照我们的输入完全的按照组分离了。因为按照组分离了，那么我们的β1就不是y2-β0，所以直接进行线性回归得到的结果是每一组自己的平均值。

这种情况下，直接进行线性回归是得不到组间对比的。我们必须引入对比矩阵作为辅助工具。

同时假如说上一种情况中，如果x = factor(c(1,2,1,2,3,3)) model.matrix(~x)时，这种时候直接进行线性回归，根据我们的公式得到的是组1和组2之间的差异，

以及组1和组3之间的差异，得不到组2和组3之间的差异。所以这种方式本身就是有缺陷的。最好就是老老实实用对比矩阵，这样一套代码的复用率很高~

最后我还有两点要说的，网上有一种说法是说类似于~0+factor()这种写法，是对0进行回归，factor是附加值，我个人不同意这种说法，举一个例子，比如说我们~a+b的含义是有两个自变量，也就是线性方程是y=β0+β1x1+β2x2，不是说不对b回归啊，同理，上面也不是说不对factor回归。

第二，我思考了一下为什么要特殊处理+0这种情况，对了大家可以试一下，带常数的只能+1和+0 如果+1是等于什么都不加，所以我大胆假设这个0和1是说要不要常数项，默认是要的，用了0就是不要，如果没有常数项，传入的又是因子时，就变成分组了。然后为什么因子会被分组呢，大家可以把因子不同的levels理解成不同的自变量，因子的levels=2说明又有两个自变量，不要把因子理解成一个多维自变量，不然会很别扭，因为没有常数项又还要保留多个自变量，当然是按照自变量分组了。而如果不带0，为什么分组会少一组是因为没必要，多一组就是多考虑一个未知数，这样计算量会变大，而实际上第一列就已经考虑了，而且还带上了常数一起算，这样计算过程能少算一部分，算是优化吧。