Variable Selection in Mixed‐Effects Location‐Scale and Location‐Shift Models¶

作者: Moritz Berger, Maria Iannario
来源: Statistics in Medicine
主题: 高维统计 / 随机矩阵
相关性: 4/10
机构绿灯: Heidelberg University（US News 前 50，免分进入精读）
链接: https://doi.org/10.1002/sim.70553

一、领域脉络与小综述¶

这个方向是什么：本子方向致力于解决“具有复杂方差结构的聚类有序响应数据”的统计建模与变量选择问题。其根本问题是如何同时刻画响应的均值水平、方差异质性（尺度）以及由聚类结构带来的随机效应，并在此基础上，利用正则化方法从大量候选变量中筛选出对上述各成分有显著影响的变量，以获得一个简约、可解释的模型。
发展脉络（history）：
- 奠基工作：标准的广义线性混合模型（GLMM）通过引入聚类水平的随机截距，处理了有序响应数据中的组内相关性。但其隐含假设“随机效应仅影响均值位置，而方差（尺度）随协变量不变”过于严格。
- 主要进展（位置-尺度模型的引入）：以 Hedeker, Mermelstein & Demirtas (2008) 和 Hedeker, Demirtas & Mermelstein (2009) 为代表，他们提出了混合效应位置-尺度模型。该模型的核心创新是允许随机效应的方差随个体或聚类水平协变量变化，从而直接刻画异方差性。本文将其扩展到有序响应（累积logit模型），并称之为“位置-移位模型”。这弥补了标准GLMM无法处理子群体中变异程度不同的缺陷。
- 当前Frontier（变量选择与正则化）：面对高维协变量，模型复杂度急剧上升。参考文献中提到了 Tutz & Schauberger (2015) 的工作，他们曾针对无随机效应的有序响应位置-尺度模型（即累积模型）提出过正则化方法。Berger & Iannario (2016) 则将其进一步扩展到一类无随机效应的位置-尺度模型中，并使用了特殊的LASSO类型惩罚。
- 本文的位置：本文明确将自己定位为上述正则化方法在“混合效应”框架下的自然推广。它填补了一个明显的缺口：没有一项现有工作能同时对混合效应位置-尺度模型或位置-移位模型的位置、尺度和随机效应三个组成部分进行并行的变量选择。
子线索聚类：
- 聚类1：变量选择方法 (LASSO及其变种)。这一簇工作（Tibshirani, 1996; 以及 Tutz & Schauberger, 2015 等）的核心是使用L1惩罚来实现稀疏解。本文采用adaptive fused LASSO (Tibshirani et al., 2005 中详述了 fused LASSO) 是因为它既能惩罚单个系数的绝对值，也能惩罚相邻类别系数之间的差值，这对于有序响应模型固有的“系数随类别平滑变化”的性质尤为关键。
- 聚类2：有序响应模型 (累积模型 vs. 连续比模型等)。本文专注于累积logit模型（McCullagh, 1980 奠定了其理论）。在混合效应设定下，位置-尺度模型（Hedeker et al., 2008, 2009）与位置-移位模型（本文引入的术语）的区分：前者允许随机效应的方差也受协变量影响，后者则假设随机效应方差恒定。这是两种处理异方差性的不同路径。
- 聚类3：混合效应模型的数值优化。GLMM的惩罚似然函数通常通过EM算法或直接优化（如牛顿-拉弗森法）求解。本文的算法（结合了坐标下降法与路径搜索）属于这一簇的变体。
这个方向在追问的核心问题与主流方法及瓶颈：
- 核心问题1: 当协变量维度Q和/或K较高时，如何稳定且高效地对GLMM进行变量选择，同时避免过拟合和模型不可识别？
- 主流方法: 使用L1系列惩罚（LASSO, adaptive LASSO, fused LASSO, elastic net等）对固定效应进行正则化。瓶颈在于，惩罚项与混合效应框架中的积分步骤（需要计算边缘似然）结合后，优化问题变得高度非凸且计算代价高昂。此外，对随机效应方差成分的正则化（如将某些随机效应方差收缩到0）尚不普遍且理论性质不清晰。
- 核心问题2: 如何为有序响应模型选择合理的“尺度”部分的协变量？标准模型假设异方差源自于某些测量误差，但实际中选择哪些变量进入位置部分、哪些进入尺度部分往往缺乏先验知识。瓶颈在于，位置和尺度部分的协变量可能高度相关（甚至相同），这会严重影响估计的稳定性和识别性。
- ⚠️ 作者的framing: 作者将缺口frame成“尽管已有正则化方法用于累积模型，但它们均未处理混合效应，特别是位置-尺度和位置-移位模型”。这使得本文成为“显然的下一步” —— 只需将已有的正则化工具移植到带有随机效应的更复杂设定上。他们淡化了“混合效应模型中的变量选择理论性质（一致性、Oracle性质）尚未被证明”这一更根本的难题，而是聚焦于 方法实现和实证表现。
张力：未见明显对立引用。所有被引的工作均逐步构建了更复杂的框架。一个值得注意的潜在张力是：一些简单的随机效应模型（仅随机截距）很难与复杂的尺度模型在拟合优度上区分，作者在实证中（模拟研究）也正是比较了两种模型（位置-尺度模型 vs. 位置-移位模型）。

二、最核心、最简单的例子 / 数学问题¶

第一步：把符号、模型、可观测数据交代清楚¶

响应变量：
- \(y_{ij}\) 是第 \(i\) 个聚类（如个体）的第 \(j\) 次观测的有序响应，其类别为 \(k = 1, \dots, K\)。
- 模型中有一个潜在的连续潜变量 \(y_{ij}^*\)，与观测值的关系为：\(y_{ij} = k\) 当且仅当 \(\theta_{k-1} < y_{ij}^* \le \theta_k\)，其中 \(\theta_k\) 是阈值参数（\(\theta_0=-\infty, \theta_K=\infty\)）。
可观测数据：
- 对于每个观测 \(y_{ij}\)，我们同时观测到：
  - 位置部分协变量：\(\mathbf{x}_{ij}\) （维度 \(p\)）。
  - 尺度部分协变量：\(\mathbf{u}_{ij}\) （维度 \(q\)）。
  - 随机效应设计矩阵：\(\mathbf{z}_{ij}\) （维度 \(r\)，通常只包含1，形成随机截距）。
- 总体数据集：\(\{\mathbf{y}_i, \mathbf{X}_i, \mathbf{U}_i, \mathbf{Z}_i\}_{i=1}^{N}\)，其中 \(\mathbf{y}_i\) 是第 \(i\) 个聚类中 \(n_i\) 个观测的反应向量。
模型：
- 混合效应位置-尺度模型可以写成：
  \[\eta_{ij} = \mathbf{x}_{ij}^T \boldsymbol{\beta} + \mathbf{z}_{ij}^T \mathbf{b}_i \quad \text{(位置线性预测子)}\]
  
  \[\tau_{ij} = \mathbf{u}_{ij}^T \boldsymbol{\gamma} \quad \text{(尺度线性预测子)}\]
  
  \[\text{logit}[P(y_{ij} \le k | \mathbf{b}_i, \mathbf{x}_{ij}, \mathbf{u}_{ij})] = \tau_{ij} (\theta_k - \eta_{ij})\]
- 未知参数：
  - \(\boldsymbol{\beta}\)：位置部分固定效应系数向量（\(p\) 维）。
  - \(\boldsymbol{\gamma}\)：尺度部分固定效应系数向量（\(q\) 维）。
  - \(\mathbf{b}_i\)：第 \(i\) 个聚类的随机效应向量（\(r\) 维），服从 \(\mathbf{b}_i \sim N(\mathbf{0}, \boldsymbol{\Sigma})\)。当采用“位置-移位模型”时，\(\tau_{ij} \equiv 1\) （即尺度部分无协变量）。
  - \(\theta_k\)：阈值参数（\(K-1\) 个）。注意模型中 \(\tau_{ij}\) 决定了尺度，它等价于在潜变量模型中的误差项方差（允许异方差）。
- 要估的对象：所有参数 \(\boldsymbol{\beta}, \boldsymbol{\gamma}, \boldsymbol{\Sigma}\)（以随机效应方差-协方差 \(\boldsymbol{\Sigma}\) 参数化，一般假设为对角阵）和阈值 \(\theta_k\)。

第二步：讲最小内核¶

本文的最小特例是移除随机效应（即混合效应退化为固定效应），且只考虑最简单的累积logit模型（位置-尺度模型，LLSM 的部分）。此时：

设定：\(y_{ij}\) 独立同分布（无聚类结构）。\(K=3\) （三个反应类别：低、中、高）。
模型：\(\text{logit}[P(y_{ij} \le k)] = \gamma_0 (\theta_k - \mathbf{x}_{ij}^T \boldsymbol{\beta})\)。为了简化，设 \(\gamma_0=1\) （位置-移位模型的特例）。
- 这就是一个标准的累积logit模型（比例优势模型）。
核心数学困难：该似然函数是凸的，且梯度、Hessian具有封闭形式，因此可以轻松估计。这不是难产的来源。
真正的困难在于选择性正则化：作者在位置部分和尺度部分都加入了变量。假设 \(p = q = 3\)，参数 \(\boldsymbol{\beta} = (\beta_1, \beta_2, \beta_3)^T\)，\(\boldsymbol{\gamma} = (\gamma_1, \gamma_2, \gamma_3)^T\)。
最小内核命题：作者想同时惩罚 \(\boldsymbol{\beta}\) 和 \(\boldsymbol{\gamma}\)，使得最终模型只有少数几个变量被选入位置部分（例如只有\(x_1\)）和尺度部分（例如只有\(u_2\)）。
关键想法：引入一个惩罚项 \(P(\boldsymbol{\beta}, \boldsymbol{\gamma}) = \lambda_\beta \sum_{m=1}^3 w_{\beta,m} |\beta_m| + \lambda_\gamma \sum_{s=1}^3 w_{\gamma,s} |\gamma_s|\)。这里 \(w\) 是自适应权重（来自adaptive LASSO），用来避免对较大系数施加过重惩罚。最小化 \(-\log L + P(\cdot)\)，其中 \(L\) 是累积logit的似然。
为什么这是最小内核：它抓住了整篇论文的核心理念：在面对一个复杂（测度可加、似然非高斯）的模型时，通过一个简单、通用的L1式惩罚项来同时实现变量选择和高维正则化，而不需要理论的识别性证明。本文的所有技术（算法、模拟）都是在这个最小内核之上增加了混合效应（非凸的似然面）而展开的。

三、这篇论文做了什么（本次重心）¶

三句话：① 本文提出了一种用于混合效应位置-尺度（LLSM）和位置-移位（LSM）模型（针对有序响应数据）的自适应融合LASSO变量选择方法。② 核心工具是 惩罚似然 (PL) 加上 坐标下降与路径搜索 (COOB) 优化算法，以同时选择位置、尺度和随机效应部分的协变量。③ 主要结论是：通过模拟和真实数据（SHARE案例）验证，该系统在产生稳定、可解释的参数估计方面有效，且正则化方法优于不进行变量选择的完整模型（在估计精度和预测性能上均有提升）。
关键设定与假设：
- 模型：累积logit位置-尺度模型（公式1，含随机效应 \(\mathbf{b}_i\)）。
- 假设 (文中未明确列为假设清单，但隐含于建模选择中):
  - 1. 随机效应分布：\(\mathbf{b}_i \sim N(\mathbf{0}, \boldsymbol{\Sigma})\)，且 \(\boldsymbol{\Sigma}\) 为对角阵（假设随机效应不相关）。
  - 1. 尺度部分参数化：\(\tau_{ij} = \exp(\mathbf{u}_{ij}^T \boldsymbol{\gamma})\) 或 \(\tau_{ij} = \exp(\mathbf{u}_{ij}^T \boldsymbol{\gamma})\) 的对数尺度线性形式（他使用了对数链接来保证方差非负）。
  - 1. 正则化：非正则化模型 假设了\(\lambda=0\)，我们不处理其识别性。
- 变量选择的具体设定：
  - 对 \(\boldsymbol{\beta}\), \(\boldsymbol{\gamma}\) 和 \(\boldsymbol{\Sigma}\) 中的对角线元素施加自适应融合LASSO惩罚。融合指对相邻阈值差进行惩罚（惩罚系数在类别间的平滑变化）。
  - 自适应权重来自非惩罚估计（如MLE）的逆。
主要结果：
- 理论结果：文中没有任何定理或渐近性质证明。所有结果都是基于模拟和案例研究的实证发现。
- 模拟研究结果 (关键):
  - 比较了3种模型（位置-尺度 Full / 位置-移位 Full / 正则化位置-尺度）。模拟中的真实数据生成自位置-尺度模型，位置部分含2个变量，尺度部分含1个变量。
  - 正则化后的模型在估计偏差 (bias) 和均方根误差 (RMSE) 上显著优于未正则化的完整模型（图3-5）。例如，对于位置部分不存在的虚假变量（\(\beta=0\)），未正则化模型估计的中位数偏差高达1.5，而正则化模型收缩到了0附近。
  - 变量选择的正确性：在所有模拟场景中，正则化模型调整了调和参数 \(\lambda\) 后，正确选择了真实位置与尺度变量的比例超过90%（表2）。
  - 随机效应方差：正则化模型能正确将不存在的随机效应方差收缩到接近0（尽管未证明其能否以概率1选择到0）。
  - 预测性能：在预测新观测的训练集-测试集划分中，正则化模型的AUC（曲线下面积）显著优于未正则化的完整模型（表3）。
- 真实数据案例 (SHARE):
  - 数据：来自欧洲健康、老龄与退休调查（SHARE），用于预测老年人的生活质量（有序5类）。
  - 变量：位置和尺度部分各使用了9个变量，包括收入、健康自评、社会参与、抑郁等；随机效应部分考虑国家层面（11个国家，随机截距）。
  - 结果：正则化模型最终保留了一个非常简约的模型。位置部分只保留了健康自评和抑郁得分；尺度部分只保留了收入和社会参与（意味着这些群体在生活质量感知上的变异更大）。随机效应部分显示不同国家间（随机截距）和不同性别内（尺度部分的随机效应？文章并未细化随机效应选择的具体结果）存在混杂。这揭示了哪些人群内部差异性更大，而不仅仅是均值更低。
- 案例分析优势：作者对比了非正则化 Full 模型（很多不显著的变量导致过拟合，解答混乱）与正则化模型（简明、可解释）。这个例子说明，在变量选择时，正则化方法提供了比手动后向选择更系统、更可复现的方式。
证明路线与技术技巧（理论型论文必写，但本论文为应用型，此处无证明路线，改为方法呈现）：
- 整体路线 (算法步骤)：
  1. 初始估计：通过标准的无惩罚GLMM程序（使用扩展后的位置-尺度模型）得到一个 \(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}},\hat{\boldsymbol{\Sigma}}\) 的初始估计。
  2. 构建自适应权重：计算每个参数的adaptive weight = \(1/|\hat{\beta}_m|\) 等。
  3. 优化惩罚似然：对惩罚后的对数似然函数 \(l(\boldsymbol{\beta}, \boldsymbol{\gamma}, \boldsymbol{\Sigma}) - \lambda_\beta P_\beta(\boldsymbol{\beta}) - \lambda_\gamma P_\gamma(\boldsymbol{\gamma})\) 进行优化。这里的关键是：
    - 使用坐标下降 (CO)：依次更新每个参数，而将其它参数视为固定。
    - 使用路径搜索 (OOB)：沿着正则化参数 \(\lambda\) 的序列（从大到小）优化，使用上一个解作为当前解的起点（热启动，warm-start），这在寻找最优 \(\lambda\) 时大大提高了计算效率。
  4. 选择 \(\lambda\): 使用BIC（贝叶斯信息准则）或交叉验证选择最优的正则化强度。
- 关键技巧:
  - 高斯-正交法：计算边缘似然时，随机效应分布的积分使用Gauss-Hermite quadrature近似。在正则化框架下，这种数值积分方法并未与参数选择解耦，而是需要反复调用，因此效率是主要关注点。
  - 自适应fused LASSO：对相邻阈值 \(\theta_k\) 的惩罚确保它们在空间中平滑，这通常是在有序响应模型中使用fused LASSO的主要动机。
真实例子与应用：已详细讲述（SHARE案例）。它想说明的核心信息是：正则化能生成一个简洁、理论通透的真实世界模型，其预测性能不差于复杂模型，且您可以根据模型简洁地讲述一个科学的叙事（哪些因素改变均值，哪些改变方差）。
🔎结论是否比证明窄：非常明显，结论比证明窄。
- 论文证明（展示）的所有内容都是基于算法+在某一个数据集上的表现。它在数值上展示了“这套方法能选对变量”，但 从未证明 它能以概率1正确选择重要变量（变量选择的一致性），也未证明其估计的渐近正态性或Oracle性质。
- 文中使用了“stable parameter estimates” (稳定的参数估计) 等语，它仅仅是基于模拟中的RMSE表现，而不是指统计意义上的稳定性（如标准误）。
- 一个潜在的风险是，用于初始估计的full（非正则化）模型在某些高维场景下根本不收敛（算法无法收敛）。作者未提及这种情况。

四、开放问题（点到为止，扎根具体语句）¶

一致性理论 (Consistency): 扎根基 §4 模拟研究 中，正则化方法的“正确变量选择比例”仅基于有限样本模拟。“在 \(p,q \to \infty\) 但稀疏性条件下，自适应融合LASSO \(\lambda\) 是否能渐近选择出所有真正重要变量并排除所有虚假变量？” 这仍是开放问题。
随机效应方差的正则化性质: 文章对随机效应方差 \(\boldsymbol{\Sigma}\) 的对角元素施加了L1惩罚。“是否能在理论上证明，该惩罚能保证真正不存在的随机效应被淘汰到0？或者这种惩罚对估计量的MSE有何影响？” 论文的模拟显示其能将不相关随机效应方差收缩到近0，但无理论保障。
尺度部分与位置部分的识别性: 当位置变量 \(\mathbf{x}\) 和尺度变量 \(\mathbf{u}\) 高度相关（甚至重合）时，模型参数 \(\boldsymbol{\beta}\) 和 \(\boldsymbol{\gamma}\) 可能无法唯一识别。“在没有任何正交性假设（如X和U的协方差结构）的情况下，该变量选择过程是否会产生虚假的路径依赖，即算法仅因为变量被排入某个组成部分而错误地将其选出？” 这论文未提及。
计算的效率瓶颈: 文中算法依赖Gauss-Hermite quadrature的反复调用并进行坐标下降。“对于更深的随机效应结构（多个水平、交叉随机效应）或更高维度（Q>100），该优化能否在可接受时间内计算？是否存在更高效的积分逼近方法（eg. Laplace近似、Hamiltonian Monte Carlo）？” 这论文仅处理了单个随机截距（国家），对于更复杂的非对角\(\boldsymbol{\Sigma}\)尚待验证。

Maintained by 陈星宇 · Homepage · Source on GitHub