前言

最近看到一篇文章，是关于利用RNA-seq来讨论两个物种之间直系同源基因的表达保守性的文章（矩阵为标准化后的）：《A comparative encyclopedia of DNA elements in the mouse genome》
其中有一张图讨论的是两个物种的直系同源基因表达量的保守性：

那么怎么理解这个残差值？作者在method部分给出了解释：

To assess the contribution of tissue and species to gene expression variation, we used a linear mixed model (LMM). Gene expression was modeled as a function of tissue and the species (considered as random factors). The LMM was implemented in the R package lme4 (ref).

The restricted maximum likelihood (REML) estimators for the random effects of tissue, species and residual variance were normalized by their sum to give the variance components.We selected genes whose fraction of variance explained by tissues or species was in their respective top quartile and higher than the other fraction to cluster the samples

由method里面我们可以了解到，作者利用lme4建立混合线性模型，并将tissue或species视为随机效应变量；然而method里面没有提到固定效应是哪一个，所以这里我们就默认为没有固定效应
另外一个信息是作者利用混合线性模型的残差值来作为指标，判断哪些基因是表达保守的，并且残差是经过标准化的

混合线性模型

作者这里是利用lme4来构造混合线性模型的，其中建模可参考：传送门 第六页有一定的介绍
其基本语法是：

mod= lmer(data = , formula = y ~ Fixed_Factor + (Random_intercept + Random_Slope | Random_Factor))

1.可以简单的把 "（）"内的内容 看为随机效应部分 ，"（）"前面的 是固定效应部分 。
2.固定效应部分不写或者写1表示只有固定效应的截距，而没有固定效应的斜率，如果赋予变量 Fixed_Factor ，则代表同时考虑固定效应的斜率和截距；
3.随机效应部分 "竖线 | "左边表示的是随机效应的斜率 ，左边写 1 即表示没有随机效应的斜率而只考虑随机效应的截距， "竖线 | "右边表示随机效应的随机因子，随机因子代表的在随机效应下的更进一步分组，例如随机效应为两组（A和B），随机因子也为两组（C和T），那么 (Random_Slope | Random_Factor) 这种写法分组就为（A-C，A-T，B-C，B-T），也就是说随机效应比较B和A的差异需要指明条件，即在 C 条件下的B和A的差异或者是T 条件下的B和A的差异，如果没有随机因子，那么随机效应系数计算就仅仅考虑B和A的差异（选择的数据包括C+T 条件下B与C+T 条件A的差异），而不加以区分C 条件还是T 条件B和A的差异（详细看后面的例子）；Random_intercept 通常写作成 1 来表示，(1+ Random_Slope | Random_Factor) 等价于 (Random_Slope | Random_Factor)，1 代表拥有截距

举个例子
例子1

dat = sleepstudy
dat
> dat
    Reaction Days Subject
1   249.5600    0       A
2   258.7047    1       A
3   250.8006    2       A
4   321.4398    3       A
5   356.8519    4       A
6   414.6901    5       A
7   382.2038    6       A
8   290.1486    7       A
9   430.5853    8       A
10  466.3535    9       A


dat$Subject = as.factor(c(rep('A',60),rep('B',60),rep('C',60)))
fm1 <- lmer(Reaction ~ Days + (1 + Subject | Days), dat)
summary(fm1)

# 随机效应的因子水平
contrasts(dat$Days)
  1 2 3 4 5 6 7 8 9
0 0 0 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 0 0
2 0 1 0 0 0 0 0 0 0
3 0 0 1 0 0 0 0 0 0
4 0 0 0 1 0 0 0 0 0
5 0 0 0 0 1 0 0 0 0
6 0 0 0 0 0 1 0 0 0
7 0 0 0 0 0 0 1 0 0
8 0 0 0 0 0 0 0 1 0
9 0 0 0 0 0 0 0 0 1

coef(fm1)
$Days
       SubjectB      SubjectC (Intercept)     Days
0  2.399840e-08  2.920663e-08    251.4051 10.46729
1  4.067253e-08  4.970416e-08    251.4051 10.46729
2 -8.451885e-09 -1.019267e-08    251.4051 10.46729
3  9.618661e-09  1.181750e-08    251.4051 10.46729
4 -9.571652e-09 -1.162255e-08    251.4051 10.46729
5  4.944233e-08  6.047095e-08    251.4051 10.46729
6 -6.091995e-09 -7.603941e-09    251.4051 10.46729
7  2.267801e-08  2.766016e-08    251.4051 10.46729
8  5.934868e-08  7.253005e-08    251.4051 10.46729
9  7.305345e-08  8.927133e-08    251.4051 10.46729

1. 例子中随机因子为Days，因此随机效应系数显示为Days，也就是说随机效应系数的求解如下，如SubjectB vs SubjectA 的系数求解为 0 天 （或者 1-9 天分别计算）的 SubjectB vs SubjectA 计算回归系数；若没有随机因子，求解SubjectB vs SubjectA则将 0-9 天的混合起来，取全部的SubjectB和SubjectA，计算SubjectB vs SubjectA 计算回归系数
2. 这里的设计矩阵（0-1）矩阵根据用户自定义，公式为：

3. 由于Subject是因子变量，这里的对照组是Days₀，因此存在SubjectB时，SubjectB=1，否则SubjectB=0，斜率2.399840e-08代表的是 Days₀ SubjectB 的 Reaction 与 Days₀ SubjectA 的 Reaction 相比较下的斜率

4. Days₀ SubjectB 的 Reaction 数据代表同时包含SubjectB和Days₀标签的数据都算；Days₀ SubjectA 的 Reaction 数据代表同时包含SubjectA和Days₀标签的数据都算；Days₂ SubjectB 的 Reaction 数据代表同时包含SubjectB和Days₂标签的数据都算；Days₂ SubjectA 的 Reaction 数据代表同时包含SubjectA和Days₂标签的数据都算
5. 由于随机效应写法为 (1 + Subject | Days)，因此

0 天的系数 SubjectB = 2.399840e-08，代表 0 天所有带有SubjectB标签的数据与0 天所有带有SubjectA标签的数据做比较所得的系数

例子2

dat = sleepstudy
dat = subset(dat,dat$Days=='0' | dat$Days=='1' | dat$Days=='2')
>dat
   Reaction Days Subject
1   249.5600    0       A
2   258.7047    1       A
3   250.8006    2       A
11  222.7339    0       A
12  205.2658    1       A
13  202.9778    2       A
21  199.0539    0       A

dat$Subject = as.factor(c(rep('A',18),rep('B',18),rep('C',18)))
fm1 <- lmer(Reaction ~ Days + (1 + Subject | Days), dat)
summary(fm1)
coef(fm1)

# 随机效应的因子水平
contrasts(dat$Subject)
 B C
A 0 0
B 1 0
C 0 1

coef(fm1)

  SubjectB  SubjectC (Intercept)     Days
0 0.1044361 0.1593605    257.5744 4.171244
1 3.2736827 4.9953630    255.5472 4.171244
2 1.9214133 2.9319144    256.4121 4.171244

1. 例子中随机因子为Days，因此随机效应系数显示为Days，也就是说随机效应系数的求解如下，如SubjectB vs SubjectA 的系数求解为 0 天 （或者 1-2 天分别计算）的 SubjectB vs SubjectA 计算回归系数；若没有随机因子，求解SubjectB vs SubjectA则将 0-2 天的混合起来，取全部的SubjectB和SubjectA，计算SubjectB vs SubjectA 计算回归系数，公式为：

2. 由于Subject是因子变量，这里的对照组是Days₀，因此存在SubjectB时，SubjectB=1，否则SubjectB=0，斜率0.1044361代表的是 Days₀ SubjectB 的 Reaction 与 Days₀ SubjectA 的 Reaction 相比较下的斜率

3. Days₀ SubjectB 的 Reaction 数据代表同时包含SubjectB和Days₀标签的数据都算；Days₀ SubjectA 的 Reaction 数据代表同时包含SubjectA和Days₀标签的数据都算；Days₂ SubjectB 的 Reaction 数据代表同时包含SubjectB和Days₂标签的数据都算；Days₂ SubjectA 的 Reaction 数据代表同时包含SubjectA和Days₂标签的数据都算
4. 由于随机效应写法为 (1 + Subject | Days)，因此

0 天的系数 SubjectB = 0.1044361，代表 0 天所有带有 SubjectB 标签的数据与 0 天所有带有 SubjectA 标签的数据做比较所得的系数

文章的例子

那么我们回顾到文章中，理想情况如下：

ddd = data.frame(gene_1 = c(10,11,22,25,32,36,
                            84,86,77,78,99,103),
                 gene_2 = c(10,11,39,40,66,68,
                            13,15,36,39,71,73),
                 species = c('h','h','h','h','h','h','m','m','m','m','m','m'),
                 tissue = c('a','a','b','b','c','c','a','a','b','b','c','c'))

## a,b,c分别表示不同的器官
## h,m 表示不同的物种

library(lme4)

fm1 <- lmer(gene_1 ~ 1 + ( 1 | species) , ddd)
coef(fm1)

$species
  (Intercept)
h    22.96185  ### 代表物种 h 的基因平均表达量
m    87.53815  ### 代表物种 m 的基因平均表达量

fm2 <- lmer(gene_1 ~ 1 + ( 1 | tissue) , ddd)

fm3 <- lmer(gene_2 ~ 1 + ( 1 | species) , ddd)

fm4 <- lmer(gene_2 ~ 1 + ( 1 | tissue) , ddd)

关于模型，我按照species 和 tissue 分别进行建模，并且仅考虑随机效率的截距，而不考虑随机效应的斜率，相当于随机效应的斜率为0。然后分别提取模型的模型随机效应的残差

gene_1 反应了该基因表达在不同物种中表达是不一样的；而gene_2反应的是该基因在不同器官里面的表达是不一样的

fm1，fm2，fm3，fm4的残差分别是115.4，1263，602.1， 5.861
对于gene_1来说，我们计算残差的百分比（即以species 和 tissue 建模的残差，除以它们的和；并且 species 为纵坐标，tissue 为横坐标），最后的结果为（0.91，0.084）；而gene_2为（0.0096，0.99）
那么对照文献里面的图

如何理解残差与保守性的关系

我们这里是利用混合线性模型进行建模，那么残差是描述模型拟合好坏的指标，残差越小，模型拟合越好；残差越大，模型拟合越差

那么如果我们仅仅以随机效应的截距建模（不考虑随机效应的斜率）：
如果我们以 species 建模，残差很小，说明模型用一条水平线（截距）去拟合模型的效果很好，也就说明该基因的表达量在两个物种中的表达“在一条水平线上”，即在不同物种中该基因表达是相同的；以 species 建模，残差很大，说明模型用一条水平线（截距）去拟合模型的效果很差，也就说明该基因的表达量在两个物种中相差很大，即在不同物种中该基因表达在不同物种中表达是不一致。