全文链接：https://tecdat.cn/?p=38816
原文出处：拓端数据部落公众号
在数据分析领域，当数据呈现出层次结构时，传统的一般线性模型（GLM）可能无法充分捕捉数据的特征。混合效应回归作为GLM的扩展，能够有效处理这类具有层次结构的数据，如聚类数据、重复测量数据和纵向数据等。本文将深入探讨混合效应回归的基本原理、关键概念、不同模型类型的差异，以及如何使用Python进行建模和分析。
混合效应回归基础
（一）定义与模型公式
混合效应回归是对一般线性模型的扩展，它考虑了数据的层次结构 。一般线性回归方程为：
其中，XX 是一个 N×pN×p 的设计矩阵，包含每个个体（NN）对于模型中每个自变量（pp）的观测值；ββ 是一个 p×1p×1 的列向量，包含模型中每个自变量的回归系数；ϵϵ 是一个 N×1N×1 的列向量，包含模型的误差（残差）。
而混合效应模型方程为：
其中，ZZ 是一个 N×qN×q 的设计矩阵，包含每个个体（NN）对于随机效应的每个协变量（qq）的观测值；uu 是一个 q×1q×1 的向量，包含矩阵 ZZ 中 qq 个协变量的随机效应。
（二）层次结构
在混合效应模型中，数据的层次结构通常用“层级”或“聚类”来描述。例如，在研究学生标准化考试成绩时，假设没有统一的课程和指导方针，数据是从不同学区的不同学校随机抽取的，每个数据行代表一个学生。这里，第1层是分析的基本单位，即学生；第2层是学校，将第1层的所有学生聚类到不同学校；第3层是学区，将第2层的学校进一步聚类。第2层及以上层级是被建模的随机效应。
图1展示了这种层次结构：
图1 数据层次结构可视化
（三）固定因素与随机因素
固定效应参数描述了整个总体中协变量与因变量之间的关系，而随机效应则特定于总体中的主体聚类。固定因素是研究感兴趣的自变量，如治疗类别、性别等；随机因素是分析单位所属的分类变量，通常定义了第2层、第3层或更高层级。例如，在上述学校研究中，学校（第2层）是随机因素，因为它是学生（第1层）的聚类变量 。变量被定义为固定效应还是随机效应，取决于研究目标和分析方法。
（四）混合效应模型类型差异
混合效应模型主要有随机截距模型、随机斜率模型和随机截距与斜率模型。随机截距模型允许基于聚类变量有不同的截距；随机斜率模型允许基于某个变量有不同的斜率；随机截距与斜率模型则同时允许基于聚类变量有不同的截距和基于某个变量有不同的斜率。
如图2所示，展示了随机截距模型和随机截距与斜率模型的差异：
图2 随机截距模型和随机截距与斜率模型差异
混合效应回归的假设与检验
（一）假设条件
- 误差独立性：各观测值的误差之间相互独立。
- 误差方差齐性：不同观测值的误差方差相等。
- 误差正态性：误差服从正态分布。
（二）假设检验
- 固定效应假设检验：
- 多个固定效应检验：
- 单个固定效应检验：
- 协方差参数似然比检验：假设嵌套模型和参考模型具有相同的固定效应，但协方差参数不同。计算参考模型和嵌套模型的 -2 REML对数似然的正差值，然后根据适当的 χ2χ2 分布查找 pp 值。
当计算的检验统计量小于指定显著 pp 值的临界值时，拒绝原假设。
Python实现混合效应回归
（一）数据准备
本研究使用的数据集，旨在比较不同窝中大鼠幼崽的出生体重。在这些窝中，一些雌性母鼠接受了“高剂量”“低剂量”和“对照”的实验处理。原设计中每个处理条件分配了10只雌性大鼠，但“高剂量”处理中有3只死亡，导致研究设计不平衡。在本研究中，窝将作为聚类变量（第2层），大鼠幼崽作为分析单位（第1层）。
- import pandas as pd
- import researchpy as rp
（二）数据探索
查看数据集中的变量信息：
分析大鼠幼崽体重基于性别和处理组的情况：
可视化体重按处理组和性别的分布：
boxplot = data.boxplot
图3 体重按处理组和性别的箱线图
（三）模型构建
- 随机截距模型
- import statsmodels.formula.api as smf
- model.summary()
移除不显著的交互项后重新建模：
"weight ~ litsize + C(treatment) + C(sex, Treatment('Male'))", data, 
计算组内相关系数（ICC）：
这表明同一窝中体重之间存在中等程度的相关性。
2. 随机斜率模型
- 随机截距和斜率不相关
- "weight ~ litsize + C(treatment) + C(sex)", data, groups= "litter",
- vc_formula = {"sex" : "0 + C(sex)"}).fit()
- 随机截距和斜率相关
- "weight ~ litsize + C(treatment) + C(sex)", data, groups= "litter",
- re_formula = "1 + C(sex)").fit()
计算随机截距和随机斜率之间的估计相关系数：
这表明体重较高的窝中，雄性大鼠幼崽往往体重也较高。
（四）假设检验
- 正态性检验
- 可视化残差的核密度估计图和Q-Q图：
图4 模型残差的KDE图
图5 模型残差的Q-Q图
- 正式的Shapiro-Wilk正态性检验：
结果显示残差的正态性假设被违反。
2. 方差齐性检验
- 可视化残差与拟合值的散点图（RVF图）和残差按窝的箱线图：
- fig = plt.figure(figsize = (16, 9))
- ax = sns.scatterplot(y = model.resid,
图6 RVF图
图7 残差按窝的箱线图
- 正式的White’s拉格朗日乘数异方差检验：
- t p-value"]
- for key, val in dict(zip(labels, het_white_res)).items():
- print(key, val)
正式检验表明方差齐性假设被违反。
结论
本文全面介绍了混合效应回归模型，从理论基础到Python实现，包括模型的构建、假设检验以及结果分析。通过对大鼠幼崽体重数据的分析，展示了混合效应回归在处理具有层次结构数据时的有效性。在实际应用中，需根据数据特点和研究目的选择合适的混合效应模型类型，并严格检验模型假设，以确保分析结果的可靠性。未来的研究可以进一步探索如何更好地处理假设违反的情况，以及将混合效应回归应用于更复杂的数据场景。
