我试图用简单的方式解释生物信息学和计算生物学中使用的一些神秘的分析技术。线性混合模型（也称线性混合效应模型）在生命科学中应用广泛，有许多教程介绍了如何在 R 中运行该模型，但有时并不清楚在概率最大化过程中究竟如何优化随机效应参数。在上一篇文章线性混合模型中，我介绍了该模型的概念，而在本教程中，我们将应用最大似然方法，从头开始推导并编代线性混合模型（LMM），也就是说，我们将使用纯R语言编码LMM，并将输出结果与lmer和lme R函数的输出结果进行比较。
让我们生成一个测试数据，结构简单为了用R代码计算简单快速，但仍然保留了线性混合建模（LMM）典型设置的所有必要元素。假设我们只有 4 个样本： 其中 2 个来自 1 号个体，另外 2 个来自 2 号个体。
此外，这 4 个点分布在两种情况下： 未治疗和治疗。假设我们测量了每个人对治疗的反应（Resp），并希望了解治疗是否导致研究中的个人产生显著反应。换句话说，我们的目标是实施与配对 t 检验类似的方法，并评估治疗的显著性。稍后，我们将把LMM 和配对 t 检验的输出结果对比起来，并证明它们实际上结果是相同。在测试数据集中，治疗列中的 0 表示 “未治疗”，1 表示 “已治疗”。首先，我们将使用不考虑数据点之间相关性的普通最小二乘法（OLS）线性回归。
library("ggplot2")
df <- data.frame(Treat = c(0, 1, 0, 1),
Resp = c(10, 25, 3, 6),
Ind = c(1, 1, 2, 2))
ggplot(df, aes(x = Treat, y = Resp, color = factor(Ind), group = Ind)) +
geom_point(size = 3) + geom_line() + labs(color = "Individual")
summary(lm(Resp ~ Treat, data = df))
虽然上述回归技术技术上是可以强行执行的，但这并不是一个好的拟合，我们在这里了一个问题。普通最小二乘法（OLS）线性回归假定所有观测点（图上的数据点）都是独立的，从而产生不相关的结果，最终产生正态分布残差。然而，我们知道图上的数据点属于 2 个个体，即每个个体有 2 个数据点。
原则上，我们可以为每个个体单独拟合一个线性模型。然而，这也不是一个好的拟合方法。我们每个人只有两个点，因此数量太少，无法对每个人进行合理的拟合。此外，正如我们之前文章所看到的，个体拟合并不能说明群体概况，因为其中一些个体的行为可能与其他个体拟合的行为相反。
ggplot(df, aes(x = Treat, y = Resp)) + geom_point() + geom_smooth(method = "lm")
ggplot(df, aes(x = Treat, y = Resp)) +
geom_smooth(method = "lm", level = 0.95) + geom_point() +
facet_wrap(~Ind, nrow = 1, ncol = 2)
与此相反，如果我们想拟合所有四个数据点，就需要以某种方式考虑到它们并非独立的事实，即其中两个数据点属于 1 号个体，两个数据点属于 2 号个体。这可以在线性混合模型（LMM）或配对检验中完成。
具有 Lmer 和 Lme 的线性混合模型
当观察结果之间不存在独立性时，我们就会使用 LMM。在我们的案例中，观察结果是在个体内部聚类的。让我们使用 LMM，对斜率和截距使用固定效应，对截距使用随机效应，这将导致在 Resp ~ Treat 公式中增加一个 (1 | Ind) 项：
library("lme4")|
summary(lmer(Resp ~ Treat + (1 | Ind), df, REML = FALSE))
这里的 REML=FALSE 仅表示我们使用的是传统的最大似然法（ML）优化，而不是限制最大似然法（我们下次再讨论 REML）。在 lmer 输出的随机效应部分，我们可以看到两个最小化参数的估计值： 残差方差对应的标准差（Std.Dev.）为 4.243，随机效应（个体间共享）方差与截距相关，标准差为 5.766。同样，在 lmer 输出的 “固定效应 ”部分，我们可以看到两个估计值： 1) 截距等于 6.5，2) 斜率/处理等于 9。因此，我们有 4 个优化参数，对应 4 个数据点。如果我们看一下 “测试数据集 ”部分的第一张图，并意识到未处理样本的两个值的平均值为 (3 + 10) / 2 =6.5，我们将其表示为 β1，而处理样本的平均值为 (6 + 25) / 2 =15.5，我们将其表示为 β2，那么固定效应的值就有意义了。后者相当于 6.5 + 9，即截距固定效应的估计值（=6.5）加上斜率固定效应的估计值（=9）。在这里，我们要注意随机效应和固定效应的精确值，因为我们将在以后推导和编码 LMM 时再现它们。
默认情况下，lme4 R包和 lmer 函数不提供 p 值等统计显著性度量，但如果您仍然希望获得 LMM 拟合的p 值，则可以使用nlme R包中的lme函数：
library("nlme")
summary(lme(Resp ~ Treat, random = ~ 1 | Ind, data=df, method="ML"))
同样，截距（StdDev = 5.766264）和残差（StdDev = 4.242649）具有随机效应，截距（值 = 6.5）和斜率/治疗（值 = 9）具有固定效应。有趣的是，固定效应的标准误差及其 t 值（Treat 的 t 值=1.5）在 lmer 和 lme 之间并不一致。然而，如果我们要求 lmer 函数中的 REML=TRUE 值，则 lme 和 lmer 的固定效应统计量（包括 t 值）相同，但随机效应统计量却不同。
summary(lmer(Resp ~ Treat + (1 | Ind), df, REML = TRUE))
这就是我们下次将介绍的最大似然 (ML) 方法和受限最大似然 (REML)方法之间的区别。
LMM 与配对 T 检验之间的关系
前面我们说过，LMM 是简单配对 t 检验的一种更复杂的形式。让我们来证明，对于我们的测试数据集，它们确实给出了相同的结果。在此过程中，我们还将了解配对 t 检验和非配对 t 检验在技术上的区别。首先，考虑到样本间的非独立性，让我们在已处理样本组和未处理样本组之间进行配对 t 检验：
t.test(df$Resp ~ df$Treat == 0, paired = TRUE)
我们可以看到，配对 t 检验报告的 t 值=1.5 和 p 值=0.3743，与使用 nlme 函数或 lmer（REML=TRUE）进行 LMM 所得结果相同。配对 t 检验统计量报告的 "差异均值 = 9"也与 lmer 和 nlme 的固定效应估计值一致，请记住，我们的处理估计值 = 9 只是处理样本与未处理样本均值之差。
那么，配对 t 检验到底是做什么的呢？那么，配对 t 检验的原理是使设置看起来像单样本 t 检验，其中一组的值被检验为显著偏离零，这是第二组均值的一种。换句话说，我们可以把配对 t 检验看成是把单个拟合的截距（见第一个图）或未处理组的平均值下移到零。简单地说，这就相当于从治疗组的 Resp 值中减去未治疗组的 Resp 值，即如下图所示将 Resp 变量转换为 Resp_std（标准化 Response），然后对 Resp_std 变量而不是 Resp 进行非配对 t 检验：
df$Resp_std[df$Treat==0] <- df$Resp[df$Treat == 0] - df$Resp[df$Treat == 0]
df$Resp_std[df$Treat==1] <- df$Resp[df$Treat == 1] - df$Resp[df$Treat == 0]
t.test(df$Resp_std ~ df$Treat == 0, paired = FALSE)
我们观察到，Treat = 0 即未治疗组的 Response 值变为 0，而治疗组（Treat=1）的 Response 值则减少了未治疗组的 Response 值。然后，我们只需使用新的 Resp_std 变量进行非配对 t 检验，其结果相当于对原始 Resp 变量进行配对 t 检验。因此，我们可以总结说，LMM 重现了配对 t 检验的结果，但允许更多的灵活性，例如，不仅可以进行两组（如 t 检验）比较，还可以进行多组比较等。
线性混合模型背后的数学
现在，让我们尝试使用我们的测试数据集推导 LMM 方程组。我们将再次查看 4 个数据点，并对处理效应、 （是固定效应）和由于两个个体内的点聚类而产生的块状结构 （实际上是随机效应贡献）进行一些数学记述。我们将用 和 参数来表示响应（Response）坐标 。
这里，β1 是个体在未处理状态下的反应，而β2 是在处理下的反应。也可以说，β1 是未治疗样本的平均值，而β2 是治疗样本的平均值。变量 u1 和 u2 是分块变量，分别代表 1 号个体和 2 号个体的特定效应。最后，ϵij ∼ N(0, σ²) 是残差误差，即我们无法建模的误差，只能尽量将其最小化，作为最大似然优化问题的目标。因此，我们可以将响应变量 写成参数β、u（即固定效应和随机效应）和 ϵ 的组合，如式（1）所示。在一般形式下，这个代数方程系统可以重写为式（2），其中索引 对应于处理效应， 描述个体效应。我们还可以用矩阵形式来表示这个方程组，即公式 (3)。因此，我们可以得出以下著名的 LMM 矩阵形式，即公式 (4)。
在这里， 被称为设计矩阵， 被称为块矩阵，它表示数据点之间的关系，即它们是否来自相关的个体，甚至像我们的情况一样来自同一个个体。值得注意的是，治疗是作为固定效应建模的，因为治疗-未治疗的水平穷尽了所有可能的治疗结果。与此相反，数据的分块结构被建模为随机效应，因为个体是从人群中抽样的，可能无法正确代表整个人群。换句话说，随机效应存在误差，即 ，而固定效应则假定没有误差。例如，性别通常被建模为固定效应，因为它通常被假定为只有两个水平（男性、女性），而生命科学中的批次效应应被建模为随机效应，因为潜在的额外实验方案或实验室会产生更多的样本间系统性差异，即许多水平，从而干扰数据分析。作为经验法则，我们可以认为固定效应不应该有很多层次，而随机效应通常是多层次的分类变量，其中的层次只代表所有可能性中的一个样本，但不是全部。
让我们继续推导数据点 的均值和方差。由于随机效应误差和残差误差均来自零均值的正态分布，而 中的非零分量来自固定效应，我们可以将 的预期值表示为等式 (5)。接下来，固定效应项的方差为零，因为固定效应被认为是无误差的，因此，对于 的方差，我们得到等式 (6)。
考虑到 和 以及 , 其中 I 是一个 4 x 4 的单位矩阵。这里，**σ² 是残差方差（建模未还原误差），σs² 是随机截距效应（跨数据点共享）方差。σs²前面的矩阵称为亲缘矩阵，由式（7）给出。亲缘矩阵对数据点之间的所有相关性进行编码。例如，有些数据点可能来自基因相关的人、地理位置接近的人、有些数据点可能来自技术复制人。这些关系都在亲缘关系矩阵中进行了编码。因此，公式 (6) 中数据点的方差-协方差矩阵的最终形式为公式 (8)。得到方差-协方差矩阵后，我们就可以继续进行需要方差-协方差的最大似然函数的优化程序。
最大似然 (ML) 原理的 LMM
我们为什么要花这么多时间推导方差-协方差矩阵，它与线性回归有什么关系？原来，拟合线性模型的整个概念，以及传统常量统计中的许多其他概念，都来自最大似然（ML）原则。为此，我们需要最大化与参数 β1, β2, σs² 和 σ² 有关的多元高斯分布函数，公式 (9)。
这里 | Σ y | 表示方差-协方差矩阵的行列式。我们看到，方差-协方差矩阵的逆矩阵和行列式明确包含在似然函数中，这就是为什么我们必须通过随机效应方差 σs ²和 残差方差σ²来计算其表达式。似然函数的最大化等同于对数似然函数的最小化，即等式 (10)。
我们需要对方差协方差矩阵的行列式、逆方差协方差矩阵以及逆方差协方差矩阵与 Y − X *β项的乘积进行繁琐的符号推导*。根据我的经验，这在 R / Python 中很难做到，但是我们可以使用Maple 进行符号计算，并推导出方差协方差矩阵的行列式和逆的表达式：
在 Maple 环境中，复杂的符号推导变得容易
使用 Maple，我们可以得到方差-协方差矩阵的行列式，如公式 (11)。接下来，公式 (10) 中对数似然的最后一项采用公式 (12) 的形式。
现在我们准备使用optim R 函数对β 1 、 β 2 、 σs ²和σ ²执行对数似然函数的数值最小化：
f<-function(x)
{
sigma <- x[1]
sigmas <- x[2]
beta1 <- x[3]
beta2 <- x[4]
y11 <- 3
y12 <- 10
y21 <- 6
y22 <- 25
-(1/2)*log(2*pi)-(1/2)*log(4*sigmas^4*sigma^4 + 4*sigmas^2*sigma^6 + sigma^8) -
(1/2)*(1/((sigma^2)*(sigma^2+2*sigmas^2)))*(((y11-beta1)^2)*(sigma^2+sigmas^2) -
2*(y11-beta1)*(y21-beta2)*(sigmas^2) +
((y21-beta2)^2)*(sigma^2+sigmas^2) +
((y12-beta1)^2)*(sigma^2+sigmas^2) -
2*(y12-beta1)*(y22-beta2)*(sigmas^2) +
((y22-beta2)^2)*(sigma^2+sigmas^2))
}
optim(par=c(1,1,1,1), f, method = "L-BFGS-B", lower=c(1,1,1,1), upper=c(10,10,10,20),
hessian = TRUE, control = list(fnscale = -1))
我们可以看到，由于我们收到“convergence = 0”消息，因此最小化算法已成功收敛。在输出中， σ =4.242640687是残差标准差，它精确地再现了 lme 和 lmer（REML = FALSE）的结果。类似地， σs =5.766281297是共享标准差，它再次精确地再现了 lme 和 lmer（REML = FALSE）函数相应的随机效应截距输出。正如预期的那样，固定效应β 1 = 6.5是未处理样本的平均值，与 lmer 和 lme 的截距固定效应估计一致。接下来， β 2 = 15.5是处理过样本的平均值，它是来自 lmer 和 lme R 函数的截距固定效应估计（= 6.5）加上斜率 / 处理固定效应估计（= 9）。
太棒了！我们通过从头开始推导和编码线性混合模型 (LMM)，成功地重现了 lmer/lme 函数的固定效应和随机效应输出！
总结
在本文中，我们学习了如何在测试数据集上推导和编码具有固定和随机效应的线性混合模型 (LMM) 。我们介绍了 LMM 与配对 t 检验之间的关系，并从lmer和lme R 函数中重现了固定和随机效应参数。
在下面的评论中，让我知道生命科学中的哪些分析技术对您来说特别神秘，我将尝试在以后的帖子中介绍它们，我们将介绍最大似然 (ML) 和 限制最大似然 (REML)方法之间的区别，敬请关注。