R语言混合线性模型、多层次模型、回归模型分析学生平均成绩GPA和可视化
原文链接：http://tecdat.cn/?p=23159
混合模型在统计学领域已经存在了很长时间。例如，标准的方差分析方法可以被看作是混合模型的特殊情况。最近，混合模型有多种应用和扩展，使其能够涵盖各种不同的数据情况。
术语
对于不熟悉的人来说，围绕混合模型的术语，特别是跨学科的术语，可能有点令人困惑。你可能遇到的关于这些类型的模型的一些术语包括。
方差分量
随机截距和斜率
随机效应
随机系数
变化的系数
截距和/或斜率作为结果
分层线性模型
多层次模型（意味着多层次的分层聚类数据
增长曲线模型（可能是Latent GCM）。
混合效应模型
都描述了混合模型的类型。混合效应，或简称混合，模型一般指固定效应和随机效应的混合。对于一般的模型，我更喜欢用 "混合模型 "或 "随机效应模型"，因为它们是更简单的术语，没有暗示具体的结构，而且后者还可以适用于许多人在使用其他术语时不会想到的扩展情况。关于混合效应，固定效应指的是人们在线性回归模型中看到的典型主效应，即混合模型的非随机部分，在某些情况下，它们被称为总体平均效应。无论如何定义，随机效应只是那些特定于观察样本的效应。本文所概述的方法主要涉及观察样本是某种分组因素水平的情况。
聚类的种类
数据可能有一个或多个聚类来源，而且聚类可能是分层的，比如聚类是嵌套在其他聚类中。一个例子是对学生进行多次的学术能力测试（重复观察嵌套在学生中，学生嵌套在学校中，学校嵌套在地区中）。在其他情况下，不存在嵌套结构。一个例子是反应时间实验，参与者执行同一组任务。虽然观察结果是在个人内部嵌套的，但观察结果也是根据任务类型进行聚类的。有人用嵌套和交叉这两个术语来区分这些情况。此外，聚类可能是平衡的或不平衡的。
在下面的内容中，我们将看到所有这些数据情况下的混合效应模型。一般来说，我们的方法将是相同的，因为这种聚类实际上更多的是数据的属性而不是模型的属性。然而，了解混合模型在处理各种数据情况时的灵活性是很重要的。
随机截距模型
下面我们展示混合模型的最简单和最常见的情况，即我们有一个单一的分组/群组结构的随机效应。这通常被称为随机截距模型。
例子：学生大学平均成绩GPA
下面我们将评估预测大学平均成绩（GPA）的因素。200名学生中的每个人都被评估了六次（前三年的每个学期），因此我们在学生中进行了观察分组。我们还有其他变量，如状态、性别和高中GPA。有些会以标签和数字的形式出现。
标准回归模型
现在说说基础模型。我们可以用几种不同的方式来展示它。首先，我们从一个标准回归开始，来确定我们的方向。
我们对截距和时间的影响有系数（b）。误差（ϵ）被假定为正态分布，平均值为0，标准差为σ。
另一种写法是强调gpa的基本数据生成过程的模型，如下。
更严格地说，GPA和μ变量有一个隐含的下标来表示每个观察值，但你也可以把它看成是单个个体在单个时间点的模型。
混合模型
描述
现在我们展示一种描述混合模型的方法，其中包括每个学生的独特效应。考虑一下下面这个单一学生的模型。这表明，学生的特定效应，即GPA的偏差只是因为该学生是谁，可以被看作是一个额外的方差来源。
我们（通常）会对学生效应做如下假设。
因此，学生效应是随机的，具体来说是正态分布，均值为零，有一定的估计标准差（τ）。换句话说，从概念上讲，这个混合模型和标准回归之间的唯一区别是学生效应，平均而言，学生效应是没有影响的，但通常会因学生的不同而有一定的变化，平均而言是τ。
如果我们重新排列，我们反而可以关注模型的系数，而不是作为额外的误差来源。
或者更简洁的说
这样一来，我们就会有针对学生的截距，因为每个人都会有自己独特的影响加到总体截距上，导致每个人的截距不同。
现在我们看到截距是正态分布的，有总体截距的平均值和一些标准差。因此，这通常被称为随机截距模型。
多层次模型
另一种显示混合模型的方式常见于多层次模型的文献中。它被更明确地显示为一个两部分的回归模型，一个在样本层面，一个在学生层面。
然而，在将第二层次的部分 "插入 "第一层次后，它与前者是相同的。
请注意，我们并没有一个针对学生情景的效应。在这种情况下，情景被说成是一个固定效应，而没有随机成分。但情况绝对不是这样的，我们后面会看到。
应用
可视化
在这里，我们绘制GPA与情景（即学期）的关系，来了解起点和趋势的变化。
plot(occasion, gpa,smooth(method = 'lm')
所有学生的路径都以淡色路径显示，10个样本以粗体显示。我们稍后要做的回归所估计的总体趋势显示为红色。有两件事很突出。一是学生在开始时有很大的变数。第二，虽然GPA的总体趋势是随着时间的推移而上升但个别学生在这个轨迹上可能会有所不同。
标准回归
因此，让我们开始吧。首先，我们来看看回归，只把时间指标作为协变量，我们把它当作数字。
lm(gpa ~ occasion)
## summary(lm)
上面的数据告诉我们，开始时，即学期为零时，平均GPA，用截距表示，是2.6。此外，随着我们从一个学期到另一个学期，我们可以预期GPA会增加大约0.11分。这将是很好的，除了我们忽略了聚类。这样做的一个副作用是，我们的标准误差是有偏差的，因此基于标准误差的统计学意义的主张会有偏差。更重要的是，我们无法探索学生效应，但学生效应很有意义。
分组回归
另一种方法是对每个学生单独进行回归。然而，这种方法有很多缺点--当有很多组的时候，它不容易被总结，通常每个组内的数据很少，无法做到这一点（如在这个案例中），而且模型是过度的背景化，意味着他们忽略了学生的共同点。我们将在后面比较这样的方法和混合模型。
运行混合模型
接下来我们运行一个混合模型，它将允许学生的特定效应。在下文中，代码看起来就像你用lm做的回归一样，但有一个额外的部分来指定组，即学生的效应。(1|student)意味着我们允许截距（用1表示）因学生而异。使用混合模型，我们可以得到与回归相同的结果，但将有更多的内容可以讨论。
library(lme4)
gpa_mixed = lmer(gpa ~ occasion + (1 | student), data = gpa)
## summary(gpa_mixed)
首先我们看到，截距和时间的系数，即在这里可以称为固定效应，与我们在标准回归中看到的相同，其解释也相同。另一方面，这里的标准误差是不同的，尽管最后我们的结论在统计意义上是相同的。请特别注意，截距的标准误差已经增加。从概念上讲，你可以认为允许每个人的随机截距允许我们获得关于个人的信息，同时认识到关于总体平均数的不确定性。
虽然我们有系数和标准误差，但你可能已经注意到，lme4并没有提供p值，这有几个原因。这有几个原因，即对于混合模型，我们基本上是在处理不同的样本量，群组内的Nc，可能因群组而异（甚至是单个观测值），以及N个总观测值，这使我们在参考分布、自由度以及如何近似 "最佳 "解决方案方面处于一种模糊的状态。其他程序自动提供p值，好像没有什么问题，而且没有告诉你他们使用哪种方法来计算p值（有几种）。此外，这些近似值在某些情况下可能非常差，或者做出可能不适合情况的假设。
然而，得到置信区间是比较直接的，如下所示7。
confint(gpa)
方差成分
与标准回归输出相比，有一点是新的，那就是学生效应的估计方差/标准差（在我们之前的公式描述中是ττ）。这告诉我们，平均而言，当我们从一个学生转移到另一个学生时，GPA会有多大的变动。换句话说，即使在根据时间点进行预测后，每个学生都有自己独特的偏差，而这个值（就标准偏差而言）是整个学生的估计平均偏差。
另一种解释方差输出的方法是注意学生方差在总数中的百分比，或0.064 / 0.122 = 52%。这也被称为类内相关，因为它也是对群内相关的估计，我们将在后面看到。
随机效应的估计
运行模型后，我们实际上可以得到学生效应的估计值。我为前五个学生展示了两种方法，既可以作为随机效应，也可以作为随机截距（即截距+随机效应）。
ef(mixed)$student %>% head(5)
coef
请注意，我们不允许学期的变化，所以它对所有的学生都是一个恒定的，也就是固定的效果。
通常情况下，我们对这些效应非常感兴趣，并希望对它们有一些不确定性的感觉。可以通过预测区间来实现。或者可以直接去看它们的图。
Interval(mixed) # 用于各种模型的预测，可能使用新的数据
sim(mixed) # 随机效应估计值的平均值、中位数和SD值
plot(mixed)) # 绘制区间估计值
下面的图是每个学生的估计随机效应及其区间估计。随机效应是正态分布，其平均值为零，由水平线表示。不包括零的区间用粗体表示。
预测
现在让我们来看看标准预测与特定群组预测的对比。与大多数R模型一样，我们可以在模型上使用预测函数。
predict(mixed, re.form=NA) %>% head()
在上面的代码中，我们指定不使用随机效应re.form=NA，因此，我们对观测值的预测和我们从标准线性模型中得到的预测差不多。
predict(mixed, re.form=NA)
predict(lm)
但是每个人都有自己独特的截距，所以让我们看看当我们加入这些信息时，预测结果有什么不同。
predict(mixed)
根据估计的学生效应，学生开始于高或低于所有学生的估计截距。下面是无条件预测与包含了前两个学生的随机截距的条件预测的直观对比。
plot(x = occasion,y = gpa, color = student,prediction, group = student,y = prediction)
我们可以看到，由于截距不同，混合模型的预测结果发生了偏移。对于这些学生来说，这种转变反映了他们相对较差的起点。
聚类层面的协变量
请注意我们把混合模型描述为一个多层次模型。
如果我们在模型中加入学生层面的协变量，例如性别，我们就会有以下结果。
其中，在插入后，我们仍然有与之前相同的模型，只是增加了一个预测因素。
因此，增加群组级协变量对我们思考模型的方式没有任何不寻常的影响。我们只是把它们添加到我们的预测变量集合中。还要注意的是，我们可以将聚类层面的协变量创建为平均值或其他一些观察层面变量的总结。当聚类代表地理单位，而观察对象是人时，这一点尤其常见。例如，我们可以将收入作为个人层面的协变量，并使用中位数来代表地理区域的整体财富。
混合模型基础知识总结
混合模型使我们能够考虑到数据中的聚类。我们更好地理解了目标变量的变异性来源。我们还得到了模型中参数的具体组别估计，使我们能够准确地了解各组之间的差异。此外，这反过来又允许我们进行特定群体的预测，从而进行更准确的预测，假设由于聚类而存在明显的变异。简而言之，即使在最简单的情况下，混合模型也有很多好处。
练习
睡眠
在这个练习中，我们将使用lme4软件包中的睡眠研究数据。以下是对它的描述。
睡眠限制研究中的受试者每天的平均反应时间。在第0天，受试者有正常的睡眠量。从那天晚上开始，他们被限制在每晚3小时的睡眠时间。观察结果代表了每天给每个受试者的一系列测试的平均反应时间（以毫秒计）。
在加载软件包后，可以按以下方式加载数据。我展示了最初的几个观察结果。
运行一个以因变量为目标变量，以天为预测变量的回归。
运行一个混合模型，用随机截距表示Subject。
解释方差分量和固定效应。
添加聚类级协变量
用GPA数据重新运行混合模型，添加聚类级协变量性别或高中GPA（highgpa），或两者。解释结果的所有方面。
在模型中加入聚类层面的协变量后，学生的方差发生了什么变化？
模拟混合模型
下面代表了一种模拟随机截距模型的简单方法。
set.seed(1234) # 复制你的结果
N = Ngroups * NperGroup
y = 2 + .5 * x + u\[groups\] + e
以上哪些代表了固定效应和随机效应？现在运行混合模型。
结果与你所期望的一致吗？
在下面的内容中，我们将改变数据的各个方面，然后在每次改变后重新运行模型，然后像以前一样进行总结并得到置信区间。对于每一项，都要特别注意结果中至少有一点变化。
首先计算类内相关系数。
.
此外，创建随机效应的密度图。
改变随机效应方差/SD和/或残差/SD，注意你对ICC的新估计，并像以前一样绘制随机效应。
将数值重置为原始值。将Ngroups改为50。你在置信区间估计中看到了什么不同？
将Ngroups重新设置为100。现在将NperGroup改为10，并再次注意到置信区间与基本条件有什么不同。