翻译学习 | 混合线性模型的思考

1. 前言

这篇文档，是为那些想了解混合线性模型的人准备的。这里面很多部分，可以在很多领域中使用。我们假定大家对一些矩阵和线性回归的理论有所了解，但是更高级的知识只有模糊的认识，希望对你有所帮助。

混合线性模型，不同的学科有不同的名称，使用不同的术语去描述同样的东西。这里试图用一种比较简明的方法进行阐述，我希望这不会让你更迷惑。

2. 介绍

混合模型是用于群集数据情况的极其有用的建模工具。拥有重复测量观测数据的数据或以其他方式聚集观测数据的数据（例如学校内的学生，地理区域内的城市）非常普遍。混合模型可以以多种方式处理此类数据，但是对于刚开始使用的术语，尤其是跨学科的术语，可能有点令人生畏。

关于混合模型，您可能会遇到一些术语：

Variance components 方差成分
Random intercepts and slopes 随机截距和斜率
Random effects 随机效应
Random coefficients 随机系数
Varying coefficients 系数变化
Intercepts and slopes-as-outcomes 拦截和结果倾斜
Hierarchical linear models 分层线性模型
Multilevel models 多层模型
Growth curve models (possibly Latent GCM) 增长曲线模型（可能是潜在的GCM）
Mixed effects models 混合效果模型

所有描述混合模型的名称，有些可能更具历史性，有些则更多地出现在特定学科中，有些则可能引用某种数据结构（例如多级群集），而另一些则是特殊情况。

混合效应或简单混合模型通常是指固定效应和随机效应的混合。我更喜欢混合模型一词，因为它很简单并且没有暗示特定的结构。

3. 标准线性模型

首先，让我们从标准线性模型开始，以熟悉该表示法。为了使事情尽可能简单，同时又可以推广到常见数据情况，我将假设一些感兴趣的变量y和一个连续/数字协变量。

yi=α+βXi+eiy_i = \alpha + \beta X_i + e_i yi=α+βXi+ei
ei∼N(0,σ2)e_i \sim \mathcal{N}(0, \sigma^2) ei∼N(0,σ2)

这里：

y 是观测值
alpha 是截距
beta是X的效应值
e 是残差，它满足平均数为0，方差为Ve的正态分布

如果写为矩阵的形式：

μ = X %*% β
y = rnorm(n, μ, σ²)

在尝试达到平衡时，我怀疑这种方法可能会在不同程度上成功或失败，但是在符号和代码之间（很多教科书演示中都缺少这种东西），我希望事情会很清楚。

4. 应用实例

让我们看一些数据，开始考虑混合模型。我将使用lme4软件包中的sleepstudy数据。以下描述来自相应的帮助文件。

睡眠剥夺研究对象每天的平均反应时间。在第0天，受试者具有正常的睡眠量。从那天晚上开始，他们每晚只能睡3个小时。观察结果代表每天对每个受试者进行的一系列测试的平均反应时间（以毫秒为单位）。

让我们使用标准的线性模型来探讨持续睡眠剥夺对反应时间的影响。

这里用线性回归模型，Days为x变量，Reaction为y变量（还有人和我一样，对因变量和自变量摸不着头脑的么，用x变量和y变量更容易理解有没有！）

# > data(sleepstudy, package='lme4')
# > slim = lm(Reaction ~ Days, data=sleepstudy)
# > summary(slim)
#
# Call:
#   lm(formula = Reaction ~ Days, data = sleepstudy)
#
# Residuals:
#   Min       1Q   Median       3Q      Max
# -110.848  -27.483    1.546   26.142  139.953
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  251.405      6.610  38.033  < 2e-16 ***
#   Days          10.467      1.238   8.454 9.89e-15 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 47.71 on 178 degrees of freedom
# Multiple R-squared:  0.2865,  Adjusted R-squared:  0.2825
# F-statistic: 71.46 on 1 and 178 DF,  p-value: 9.894e-15

在天数为正的情况下，我们看到更多的睡眠剥夺会导致反应时间增加/变慢。但是让我们绘制数据。

这告诉我们什么？黑线是我们当前模型的建议，即假设每个人的出发点和轨迹都相同。但是，我们看到对象的起点可能相差100毫秒之多。此外，虽然斜率通常为正，但有些斜率随时间变化很小甚至没有变化。换句话说，个人的截距和坡度都有明显的变化。我们将在最后讨论这些数据。

5. 所有可能的混线性模型分析这个数据

因此，我们要考虑数据的集群性质。与其像上面的SLiM中那样忽略聚类，不如考虑为每个人运行完全独立的回归。但是，这些模型通常只需要很少的数据就可以运行，并且会被过度上下文化。正如我们将看到的，混合模型将允许每个人的随机截距和斜率，并在不因个人而异的情况下考虑聚类。

如何描述这个模型？事实证明，它可以并且以多种方式显示，具体取决于您正在查看的文本或文章。以下内容受Gelman＆Hill（2007）的启发，他们展示了编写混合模型的五种方法。为简单起见，我们通常只关注随机截距模型，但有时会超出该范围。前几对公式仅需要了解标准回归模型，但其他模型描述则需要更多知识。

5.1 Mixed Model 1a: Allowing coefficients to vary across groups

在上面，c簇内的每个观测i都有一个截距α，这取决于它所属的c簇。αc假设为正态分布，平均μα和方差τ2。μα是我们在SLiM方法中看到的总截距。e是SLiM中描述的正态分布的平均零。对于每一个模型描述，我将注意到一个主要的参考，在那里人们可以看到它的形式与特定的文本或文章几乎相同。它将不是唯一一个这样做的引用，但至少它应该是一个提供一些额外视角的引用。

5.2 Mixed Model 1b: Multilevel model

5.3 Mixed Model 2: Combining separate local regressions

5.4 Mixed Model 3a: Design matrix for random component

5.5 Mixed Model 3b: Design matrix again

5.6 Mixed Model 3c: General notation

5.7 Mixed Model 4a: Regression with multiple error terms

5.8 Mixed Model 4b: Conditional vs. marginal model

5.9 Mixed Model 5b: Multivariate normal model

5.10 Mixed Model 6: Penalized regression

5.11 Mixed Model 7: Bayesian mixed model

6 模拟数据运行混合模型

这里，设置：
Va = 0.50.5 = 0.25
Ve = 11 =1
随机因子blup：gamma_
截距：3
固定因子blue：0.75

# setup
set.seed(1234)
nclus = 100                                # number of groups
clus = factor(rep(1:nclus, each=10))       # cluster variable
n = length(clus)                           # total n# parameters
sigma = 1                                  # residual sd
tau = .5                                   # re sd
gamma_ = rnorm(nclus, mean=0, sd=tau)      # random effects
e = rnorm(n, mean=0, sd=sigma)             # residual error
intercept = 3                              # fixed effects
b1 = .75# data
x = rnorm(n)                               # covariate
y = intercept + b1*x + gamma_[clus] + e    # see model 1
d = data.frame(x, y, clus=clus)
head(d)
str(d)

使用lme4包运行

library(lme4)
lmeMod = lmer(y ~ x + (1|clus), data=d)
summary(lmeMod)

估算的结果可以看出：
Va = 0.224，和0.25比较接近
Ve = 0.97，和1比较接近
blue：0.799，和0.75接近
截距：2.9，和3接近

提取blup值：

asreml代码

> library(asreml)
> mod1 = asreml(y ~ x, random = ~ clus, residual = ~ idv(units),data=d)
Model fitted using the sigma parameterization.
ASReml 4.1.0 Wed Apr  5 16:34:50 2020LogLik        Sigma2     DF     wall    cpu1     -3817.282           1.0    998 16:34:50    0.02     -2862.495           1.0    998 16:34:50    0.03     -1811.528           1.0    998 16:34:50    0.04     -1082.178           1.0    998 16:34:50    0.05      -688.386           1.0    998 16:34:50    0.06      -572.668           1.0    998 16:34:50    0.07      -553.615           1.0    998 16:34:50    0.08      -552.690           1.0    998 16:34:50    0.09      -552.687           1.0    998 16:34:50    0.0
> summary(mod1)$varcompcomponent  std.error   z.ratio bound %ch
clus        0.2247008 0.04605034  4.879461     P   0
units!units 0.9755909 0.04601438 21.201871     P   0
units!R     1.0000000         NA        NA     F   0
> coef(mod1)$fixedeffect
x           0.7994379
(Intercept) 2.9008683

结果是一致的

比较设定的blup值和计算的blup值的相关系数：

> cor(gamma_,coef(mod1)$random)effect
[1,] 0.838686

7 sleepstudy数据运行混合模型

sleepMod = lmer(Reaction ~ Days + (Days|Subject), data=sleepstudy)
summary(sleepMod)

> # sleepstudy 数据
> sleepMod = lmer(Reaction ~ Days + (Days|Subject), data=sleepstudy)
> summary(sleepMod)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)Data: sleepstudyREML criterion at convergence: 1743.6Scaled residuals: Min      1Q  Median      3Q     Max
-3.9536 -0.4634  0.0231  0.4633  5.1793 Random effects:Groups   Name        Variance Std.Dev. CorrSubject  (Intercept) 611.90   24.737       Days         35.08    5.923   0.07Residual             654.94   25.592
Number of obs: 180, groups:  Subject, 18Fixed effects:Estimate Std. Error t value
(Intercept)  251.405      6.824  36.843
Days          10.467      1.546   6.771Correlation of Fixed Effects:(Intr)
Days -0.138

查看每个subject的效应值以及截距：

> ranef(sleepMod)
$Subject(Intercept)        Days
308   2.2575329   9.1992737
309 -40.3942719  -8.6205161
310 -38.9563542  -5.4495796
330  23.6888704  -4.8141448
331  22.2585409  -3.0696766
332   9.0387625  -0.2720535
333  16.8389833  -0.2233978
334  -7.2320462   1.0745075
335  -0.3326901 -10.7524799
337  34.8865253   8.6290208
349 -25.2080191   1.1730997
350 -13.0694180   6.6142185
351   4.5777099  -3.0152825
352  20.8614523   3.5364062
369   3.2750882   0.8722876
370 -25.6110745   4.8222518
371   0.8070591  -0.9881730
372  12.3133491   1.2842380with conditional variances for “Subject”

如果我们对其进行作图：

我们可以看到混合模型的好处，因为我们会有结合了个体特定影响的预测，预测的更准确。

8 其它主题

我将简要提及其他一些主题，但这些主题不会改变到目前为止讨论的一般方法。

增加分组的协变量（Cluster level covariates ）
注意随机因子是镶嵌结构，还是交互结构
你可能注意lme4包中没有给出p-value值，软件不会直接给出（除非用的是贝叶斯框架），其它软件包给出p-value，不一定说明他就是正确的。
随机因子有关系矩阵？响应变量是二分类？还有很多问题需要考虑，有些数据不适合用混合模型去分析

9. 汇总

在模型描述和代码之间，希望您对标准的混合模型框架有了很好的了解。混合模型是对标准glm的非常灵活的扩展，可以直接与加性模型，空间模型和其他模型建立联系，因此可以将它们带到很远。我可以说在lme4，mgcv和brms之间，将有很多很多方法可以以多种方式浏览其数据。祝您研究顺利！

10. 参考文献

Bates, Douglas, Martin Mächler, Benjamin Bolker, and Steven Walker. 2015. “Fitting Linear Mixed-Effects Models Using Lme4.”
Fahrmeir, Ludwig, Thomas Kneib, Stefan Lang, and Brian Marx. 2013. “Regression.”
Gelman, Andrew, and Jennifer Hill. 2007. “Data Analysis Using Regression and Multilevel/Hierarchical Models.”
Gelman, Andrew, John Carlin, Hal Stern, David Dunson, Aki Vehtari, and Donald Rubin. 2013. “Bayesian Data Analysis.”
Wood, Simon. 2006. “Generalized Additive Models.”

英文原文：https://m-clark.github.io/docs/mixedModels/mixedModels.html#standard_linear_model

个人公众号：育种数据分析之放飞自我