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和一个连续/数字协变量。
y
i
=
α
+
β
X
i
+
e
i
y_i = \alpha + \beta X_i + e_i
yi=α+βXi+ei
e
i
∼
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 2020
LogLik Sigma2 DF wall cpu
1 -3817.282 1.0 998 16:34:50 0.0
2 -2862.495 1.0 998 16:34:50 0.0
3 -1811.528 1.0 998 16:34:50 0.0
4 -1082.178 1.0 998 16:34:50 0.0
5 -688.386 1.0 998 16:34:50 0.0
6 -572.668 1.0 998 16:34:50 0.0
7 -553.615 1.0 998 16:34:50 0.0
8 -552.690 1.0 998 16:34:50 0.0
9 -552.687 1.0 998 16:34:50 0.0
> summary(mod1)$varcomp
component 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)$fixed
effect
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: sleepstudy
REML criterion at convergence: 1743.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.9536 -0.4634 0.0231 0.4633 5.1793
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 611.90 24.737
Days 35.08 5.923 0.07
Residual 654.94 25.592
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 251.405 6.824 36.843
Days 10.467 1.546 6.771
Correlation 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.2842380
with 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
个人公众号：育种数据分析之放飞自我
