具有不完整数据的混合模型的识别与不确定性量化，在化学计量学应用中的研究

《Expert Systems with Applications》：Identification and uncertainty quantification of hybrid models with incomplete data for chemometric applications

　　**中文翻译：**

**洪日聪、王波、陈涛**
石河子大学理学院，石河子，832003，中国

**摘要：**
混合模型结合了多种建模方法，在准确捕捉化学工程中的复杂系统方面至关重要。这些模型中的参数识别和不确定性量化对其预测性能和泛化能力有显著影响。在本文中，我们提出了一种新的混合模型框架，旨在处理不完整的数据集——这是化学计量学中常见的情况。为了模型识别和不确定性量化，我们开发并研究了两种方法：蒙特卡洛期望最大化（MCEM）算法和插件方法。通过广泛的模拟以及一个真实的化学物质分配系数实验数据集，我们证明了所提出框架的有效性，并评估了这四种方法在参数估计、预测准确性和不确定性量化方面的相对性能。

**1. 引言：**
许多化学和生物化学应用中的机制模型本质上都是混合和分层的，例如参见Pearson（1894年）、Pearson（1901年）、Chen、Lian和Han（2008年）以及Wang、Chen、Lian和Han（2010年）的研究。混合模型是一种统计方法，用于描述更大群体中的多个子群体。在统计学中，高层模型（或高级模型）通常代表整个群体，其特征是质量守恒、能量守恒、动量守恒等方程，以及反应和传输机制（参见Ismail（1989年）、Chen等人（2013年）和Zendehboudi、Rezaei和Lohi（2018年）的研究）。这些方程通常包含需要直接确定或通过其他已知/可测变量（如反应物浓度、温度、分子量）的经验关系来确定的参数。低层模型代表子群体（Chen、Lian和Kattou，2016年）。

第4节详细给出了一个 motivate 示例，该示例涉及对人类皮肤角质层（SC）与水（Ksc）之间化学物质分配系数的建模（Wang等人，2010年）。分配系数定义为两种相（在这种情况下是SC和水）在界面处达到平衡时的溶质浓度比。近年来，该领域的进步体现在开发了更复杂的模型，这些模型将经验数据与详细的分子理解相结合，利用了计算化学和机器学习的技术（Xia、Chen和Zhang，2023年）。这些先进的方法提高了分配系数预测的准确性和可解释性，便于其在药物开发和环境科学等领域的应用（Bannan、Calabró、Kyu和Mobley，2016年；Mackay、Celsie和Parnis，2016年）。

由于SC主要由蛋白质（存在于死去的角质形成细胞中）、脂质和水组成，Ksc 可以表示为：
$$
K_{sc} = \phi_{\text{pro}} K_{\text{pro}} + \phi_{\text{lip}} K_{\text{lip}} + \phi_{\text{w}} K_{\text{w}} + \epsilon_{\text{sc}}
$$
其中 $ K_{sc} $ 表示化学物质的分配系数，$\phi_{\text{pro}}$、$\phi_{\text{lip}}$ 和 $\phi_{\text{w}}$ 分别表示体积分数，$\epsilon_{\text{sc}}$ 表示模型中的缺陷。体积分数可以在实验中直接测量，显然对于任何化学物质都有 $ K_{w} = 1 $。然而，$ K_{\text{pro}} $ 和 $ K_{\text{lip}} $ 的值受到化学物质疏水性的影响，必须根据实验数据与辛醇-水标准分配系数（$ K_{\text{ow}} $）相关联：
$$
\begin{align*}
K_{sc} &= \frac{\phi_{\text{pro}} K_{\text{pro}} \exp\left(-\theta_1 x + \theta_2 (y - z)\right)}{\phi_{\text{lip}} K_{\text{lip}} \exp\left(-\theta_1 y + \theta_2 z\right) + \epsilon_{\text{sc}} \\
&= \frac{K_{\text{ow}} \left(\phi_{\text{pro}} \exp\left(-\theta_1 x\right) \phi_{\text{lip}} \exp\left(-\theta_2 z\right) + \epsilon_{\text{sc}}\right)}{\phi_{\text{w}} \exp\left(-\theta_1 y\right) \exp\left(-\theta_2 z\right)}
\end{align*}
$$
其中 $\theta_1 $ 和 $\theta_2 $ 包含模型参数，$ h_{\text{pro}}(x, y, z)$ 和 $ h_{\text{lip}}(x, y, z)$ 是已知的模型函数，$\epsilon_{\text{pro}}$ 和 $\epsilon_{\text{lip}}$ 表示模型中的缺陷。

在上述实验中，我们的主要关注点是模型参数的识别及其不确定性的量化。这确保了可以做出可靠的预测，并确定相关的预测不确定性。了解预测不确定性对于许多基于模型的应用（如风险评估、稳健优化和控制）至关重要。不确定性量化已成为统计学、机器学习和工程建模中的一个重要课题，因为它能够评估模型预测的可靠性和稳健性。为此目的已经开发了多种方法，包括贝叶斯校准和复杂模型的概率推断方法（Gelman、Carlin、Stern和Rubin，1995年；Higdon、Kennedy和Cavendish、Cafeo、Ryne，2004年；Kennedy和O’Hagan，2001年），用于传播不确定性的蒙特卡洛和基于抽样的技术（Robert、Casella和Casella，1999年；Sullivan，2015年），以及高斯过程回归等概率替代建模方法（Williams和Rasmussen，2006年）。这些方法已成功应用于模型校准、带有不确定性范围的预测以及许多科学和工程应用中的风险决策。模型识别和不确定性分析在过程系统工程和化学计量学领域得到了广泛研究（参见Chen等人，2008年的综述）。最近在不确定性量化方面的进展进一步提高了复杂动态和工程系统预测模型的可靠性（Cheng、Quilodrán-Casas、Ouala、Farchi、Liu和Tandeo等人，2023年；Majda和Qi，2018年；Shi、Wei和Feng，2025年）。

然而，大多数现有研究集中在单层模型上，或者假设可以获取完整的数据集进行参数估计。在混合模型的背景下，特别是在使用多层数据集进行模型识别和不确定性量化时，这一主题的探索仍然有限。实际上，许多实际建模问题涉及混合或分层结构，其中多个模型层级是耦合的，并且某些变量的观测值只是部分的。这种情况在参数识别和不确定性传播中引入了额外的挑战。以分配系数模型为例，虽然存在单独的 Ksc、Klip 和 Kpro 的实验数据，但在Wang等人（2010年）的研究中，没有使用 Kpro 数据，因为（i）Ksc 和 Klip 数据理论上足以识别模型；（ii）Ksc 的数据比 Kpro 的数据更多，因为实验测量与SC混合物的分配比与分离出的蛋白质的分配要容易得多。相反，当高层模型的特定数据较为稀少时（例如某些生物反应器建模的情况），常见的方法是仅使用低层数据来识别低层模型，然后将低层模型插入到高层模型中。这两种方法都不理想，因为它们没有利用完整的数据集。

本文通过开发一个统一的框架来解决这一问题，该框架能够在数据不完整和多层数据存在的情况下识别混合模型，并量化参数和预测不确定性。这种方法解决了混合模型的固有复杂性，标志着与传统方法的显著区别，旨在更全面地利用数据以提高模型预测的准确性和可靠性。从 motivative 示例出发，我们将提出一个相关模型的一般问题表述，然后提出两种模型识别和不确定性量化方法。第一种方法将问题框架在最大似然估计框架内。我们将介绍期望最大化（EM）算法，特别是蒙特卡洛EM（MCEM）算法（Wei和Tanner，1990年），用于模型识别和不确定性量化。虽然EM算法广泛用于涉及潜变量或不完整数据的问题，但当底层模型结构复杂时，期望步骤可能在分析上变得难以处理。在本研究中考虑的混合模型中，一些变量不是直接观测到的，似然函数涉及多个耦合的模型组件。MCEM算法通过使用蒙特卡洛抽样来近似期望步骤，从而允许参数估计过程在处理不完整的观测值和混合模型中的复杂依赖结构时保持计算可行性。因此，MCEM为本文解决的识别问题提供了灵活且实用的解决方案。第二种更直接的方法是将低层模型直接插入到高层模型中，允许使用所有数据来估计模型参数。这种方法扩展并反映了大多数现有化学计量学中关于混合和分层模型的研究实践。然而，这种方法可能导致模型不确定性的估计失真，因为它忽略了较低层模型的缺陷（如方程（2）和（3）中的 $\epsilon_{\text{pro}}$ 和 $\epsilon_{\text{lip}}$ 所捕获的）。本研究旨在揭示此类情况所面临的挑战和可能的解决方案。尽管研究将以皮肤SC中的分配系数作为具体应用，但所提出的方法是通用的。

**本工作的主要贡献包括：**
（1）提出了一个系统框架，用于在数据集不完整且某些内部变量未直接观测时识别混合模型；
（2）开发了一种统一的不确定性量化程序，用于评估参数不确定性并将其传播到混合建模结构中的模型预测；
（3）所提出的方法允许使用异构和部分观测的数据集进行实际模型校准，这在现实世界的工程和数据驱动建模应用中很常见；
（4）建立了关于模型可识别性和插件参数估计器渐近偏差的理论结果。这些特点使所提出的框架不同于假设观测值完整或单层模型结构的传统参数估计和不确定性分析方法。

**本文的其余部分组织如下：**
第2节提出了相关模型的一般问题表述，然后提出了一种贝叶斯模型识别方法，该方法将量化模型参数的不确定性，包括MCEM和插件方法。第3节通过模拟研究评估了这些方法，并在第4节将其应用于一个真实的化学数据集。最后，第5节对本文进行了一些讨论。

**2. 问题和方法论：**
**2.1. 问题表述：**
在我们的研究中，我们使用一个通用的数学模型来表示感兴趣变量之间的关系。具体来说，我们考虑以下方程组：
$$
\begin{align*}
y \in \mathbb{R}^d & \text{表示} d \text{维的高层输出变量} \\
z \in \mathbb{R}^d & \text{表示} d \text{维的底层输出变量} \\
x \in \mathbb{R}^d & \text{表示} d \text{维的输入变量} \\
\theta & \text{表示} \text{需要估计的参数向量}
\end{align*}
$$
函数 $ f(x, y, z) $ 和 $ h(x, y, z) $ 描述了这些变量之间的关系，并假设是已知的。项 $\epsilon$ 和 $\zeta$ 捕捉了关系的随机性质，考虑了测量噪声和其他不确定性来源。假设 $\epsilon$ 和 $\zeta$ 是独立的、均值为零的正态分布随机变量。它们的方差由协方差矩阵 $\Sigma_y$ 和 $\Sigma_z$ 表示，允许异方差性，并为建模变量之间的关系相关的不确定性提供了一个灵活的框架。为了简洁起见，我们假设 $\Sigma_y$ 和 $\Sigma_z$ 是对角的，这意味着在给定模型结构的情况下残差组分是条件独立的。这个假设并不意味着高层和低层变量之间是独立的，因为依赖性通过（4）和（5）中的分层结构明确捕获。该框架可以扩展以允许完整的协方差矩阵，而不会遇到数学难题。然而，在本工作中考虑的不完整和不平衡数据设置中，引入不受限制的协方差结构会显著增加参数的数量，并可能影响可识别性和数值稳定性。这里采用的对角规范代表了灵活性和统计稳定性之间的折中。我们假设 $\theta$、$\Sigma_y$ 和 $\Sigma_z$ 是独立的，它们的先验分布分别表示为 $ p(\theta)$、$ p(\Sigma_y)$ 和 $ p(\Sigma_z)$。

**2.2. 模型识别：**
在实践中，必须考虑分析中可能同时遇到的不同类型的数据点。特别是，数据点可以是：（i）包含所有变量值的完整数据；（ii）特定于低层模型的数据；（iii）特定于高层模型的数据。我们首先为不完全观测下的混合模型建立了一个可识别性结果。

**定理1：**
考虑混合模型（4）和（5）。假设：(A1) x 的分布支持集包含 R^d 中的一个非空开集。(A2) 下层回归函数 h(x, θ) 是可辨识的，即 h(x,θ1)=h(x,θ2) 几乎必然成立意味着 θ1=θ2。(A3) 对于固定的 z，上层回归函数 f(x, z, θ) 在 θ 方上是可辨识的，即 f(x,z,θ1)=f(x,z,θ2) 几乎必然成立意味着 θ1=θ2。那么，(i) 如果数据集包含完整的观测值 (x, z, y)，参数 θ 是可辨识的；(ii) 如果数据集只包含下层观测值 (x, z)，无论是否有上层观测值，参数 θ 都是可辨识的；(iii) 如果数据集不包含下层观测值 (x, z)，参数 θ 通常是不可辨识的。证明见 A 部分。

对于给定的完整数据点，用通用符号 {x_nc, z_nc, y_nc} 表示，θ、Σ_y 和 Σ_z 的后验分布由下式给出：
$$
p(\theta, \Sigma_y, \Sigma_z | x_nc, z_nc, y_nc) \propto p(y_nc | x_nc, \theta, \Sigma_z) \cdot p(\theta) \cdot p(\Sigma_y) \cdot p(\Sigma_z)
$$
当数据点仅与下层模型相关时，它将只包含变量 {x, z}。对于这样的数据点 {x_nl, z_nl}，θ、Σ_y 和 Σ_z 的后验分布变为：
$$
p(\theta, \Sigma_y, \Sigma_z | x_nl, z_nl) \propto p(z_nl | x_nl, \theta, \Sigma_z) \cdot p(\theta) \cdot p(\Sigma_y) \cdot p(\Sigma_z)
$$
另一方面，如果数据点仅与上层模型相关，则它将只包含变量 {x, y}。在这种情况下，变量 z 未被观测到，我们必须将其视为缺失变量。为了解决这个问题，我们在下一小节中提出了两种方法。

设 N_c、N_l 和 N_u 分别为完整数据、下层特定数据和上层特定数据的数据点数量，用 D 表示所有可用的数据集合。那么 θ、Σ_y 和 Σ_z 的后验分布由下式给出：
$$
p(\theta, \Sigma_y, \Sigma_z | D) \propto \prod_{i=1}^{N_c} p(y_nc | x_nc, z_nc, \theta, \Sigma_y) \cdot \prod_{i=1}^{N_l} p(z_nl | x_nl, \theta, \Sigma_z) \cdot \prod_{i=1}^{N_u} p(y_nu | x_nu, \theta, \Sigma_y, \Sigma_z)
$$
虽然所提出的框架可以应用于任何先验分布，但在本研究中，我们假设所有 θ、Σ_y 和 Σ_z 的先验都是非信息性的，以减少主观先验假设的影响，并允许推断过程主要由观测数据驱动。当关于参数的可靠先验知识不可用时，这种选择特别合适。在非信息性先验下，后验分布主要由观测数据集中的似然函数决定。因此，参数估计反映了数据中包含的信息，同时避免了来自强先验假设的额外偏差。

然后方程 (8) 可简化为：
$$
p(\theta, \Sigma_y, \Sigma_z | D) \propto \prod_{i=1}^{N_c} L_c \cdot \prod_{i=1}^{N_l} L_u
$$
其中 L_c、L_l 和 L_u 分别是完整数据、下层特定数据和上层特定数据的似然函数，它们由下式给出：
$$
L_c = \prod_{i=1}^{N_c} p(y_nc | x_nc, z_nc, \theta, \Sigma_y) \cdot p(z_nc | x_nc, \theta, \Sigma_z)
$$
$$
L_l = \prod_{i=1}^{N_l} p(z_nl | x_nl, \theta, \Sigma_z)
$$
$$
L_u = \prod_{i=1}^{N_u} p(y_nu | x_nu, \theta, \Sigma_y, \Sigma_z)
$$

为了估计参数 θ、Σ_y 和 Σ_z，我们使用最大后验概率 (MAP) 方法，即最大化 L_f 关于 θ、Σ_y 和 Σ_z。然而，由于正态性假设，L_c 和 L_l 容易计算，而 L_u 涉及到潜在变量 z_n，因此不容易获得。我们考虑两种方法来解决这个问题。

2.2.1 蒙特卡洛 EM 方法
给定一个上层特定数据点 {x_nu, y_nu}，我们将相应的 z_nu 视为潜在变量，并采用 EM 算法进行参数估计。如果 z_nu 是已知的，完整数据的对数似然将是：
$$
L_c = \sum_{i=1}^{N_u} p(z_nu, y_nu | x_nu, \theta, \Sigma_y, \Sigma_z)
$$
EM 算法中的期望步骤需要计算 〈z_nu | x_nu, y_nu, θ, Σ_y, \Sigma_z〉，由于底层模型的非线性，这没有封闭形式。这里，?·? 表示期望操作。在这种情况下，〈z_nu | x_nu, y_nu, θ, Σ_y, Σ_z〉 表示在 x_nu, y_nu, θ, Σ_y, Σ_z 给定的条件下 z_nu 的条件期望。为了解决这个问题，我们采用了一种称为 MCEM 算法的 EM 算法的变体，该算法利用蒙特卡洛抽样技术来近似所需的期望（Gelman, Carlin, Stern, Rubin, 1995; Robert, Casella, Casella, 1999）。

在 MCEM 算法中，我们使用重要性抽样技术，该技术为模拟中的不同结果分配权重，其中 p(z_nu | x_nu, θ, Σ_z) 作为重要性密度。对于第 i 个样本，记为 z_ni，我们分配以下重要性权重：
$$
w_ni = \frac{p(z_ni, y_nu, θ, \Sigma_y, \Sigma_z | x_nu)}{p(z_ni | x_nu, θ, \Sigma_z)}
$$
抽样后，需要将权重规范化以确保它们总和为 1。设 M 为抽样过程中生成的总样本数，则：
$$
w_ni = \frac{w_ni}{\sum_{i=1}^{M} w_ni}
$$
如果需要，根据它们的权重对样本进行有放回抽样可以得到权重相等的样本（1/M）。

因此，E 步骤为：
$$
\langle L_c \rangle = \sum_{i=1}^{N_u} \sum_{i=1}^{M} w_ni \cdot \ln(p(y_nu | x_nu, z_ni, \theta, \Sigma_y) + \ln(p(z_ni | x_nu, θ, \Sigma_z))
$$
在 M 步骤中，我们使用梯度下降算法等来最大化以下表达式关于 θ、Σ_y 和 Σ_z：
$$
L_f = L_c + L_l + \langle L_c \rangle
$$
其中 L_c = \ln L_c 和 L_l = \ln L_l。关于 θ 的梯度分别为：
$$
\frac{\partial L_c}{\partial \theta} = -\sum_{i=1}^{N_c} \left[ \ync - f(x_nc, z_nc, \theta) \right]^T \Sigma_y^{-1} \cdot \frac{\partial f(x_nc, z_nc, \theta)}{\partial \theta} + \left[ z_nc - h(x_nc, \theta) \right]^T \Sigma_z^{-1} \cdot \frac{\partial h(x_nc, \theta)}{\partial \theta}
$$
$$
\frac{\partial L_l}{\partial \theta} = -\sum_{i=1}^{N_l} \left[ z_nl - h(x_nl, \theta) \right]^T \Sigma_z^{-1} \cdot \frac{\partial h(x_nl, \theta)}{\partial \theta}
$$
$$
\frac{\partial \langle L_c \rangle}{\partial \theta} = -\sum_{i=1}^{N_u} \sum_{i=1}^{M} w_ni \cdot \left[ y_nu - f(x_nu, z_ni, \theta) \right]^T \Sigma_y^{-1} \cdot \frac{\partial f(x_nu, z_ni, \theta)}{\partial \theta} + \left[ z_ni - h(x_nu, \theta) \right]^T \Sigma_z^{-1} \cdot \frac{\partial h(x_nu, \theta)}{\partial \theta}
$$
如果我们定义 Σ_z^{-1} = \diag(\beta_1, \beta_2, \ldots, \beta_d) 和 Σ_y^{-1} = \diag(\alpha_1, \alpha_2, \ldots, \alpha_d)，则 〈L_c〉 关于方差的梯度分别为：
$$
\frac{\partial \langle L_c \rangle}{\partial \beta_d} = \frac{1}{2} \sum_{i=1}^{N_u} \sum_{i=1}^{M} w_ni \cdot \left[ -z_nd_i - h(x_nu, \theta) \right]^2
$$
$$
\frac{\partial \langle L_c \rangle}{\partial \alpha_d} = \frac{1}{2} \sum_{i=1}^{N_u} \sum_{i=1}^{M} w_ni \cdot \left[ -y_ndu - f(x_nu, z_ni, \theta) \right]^2
$$
在这里，下标 ‘nd’ 表示第 n 个样本的第 d 个分量，hd(·, ·, ·) 和 fd(·, ·, ·) 分别是 h(·, ·, ·) 和 f(·, ·, ·) 的第 d 个分量。L_c 和 L_l 关于方差的梯度可以类似地推导出来。注意 L_l 不涉及 Σ_y。实际上，如果我们令 ?〈L_c〉?β_d=0 和 ?〈L_c〉?α_d=0，则 β_d 和 α_d 可以通过封闭形式解得到。因此，在实践中，M 步骤可以通过以下步骤执行：(1) 对于固定的 {β_d}d=1dz 和 {α_d}d=1dy，最大化 (14) 关于 θ；(2) 给定 θ，计算 {β_d}d=1dz 和 {α_d}d=1dy。

MCEM 算法的收敛性质在文献中得到了广泛研究（Wei & Tanner, 1990）。在实际实现中，期望步骤使用蒙特卡洛抽样来近似，并且参数估计值会迭代更新，直到满足收敛标准。在本研究中，蒙特卡洛样本量选择得足够大，并通过检查连续迭代中参数估计值和相应对数似然值的稳定性来监控收敛性。当参数估计值或对数似然的相对变化低于预定义的阈值时，算法终止。这种实际的收敛性评估确保了参数估计的可靠性，同时保持了计算的可行性。

2.2.2 插值方法
作为 MCEM 方法的替代方法，另一种常用于皮肤 SC 中的分区系数建模的方法是插值方法。这种方法更直接，因为它不需要抽样或近似，因此计算速度更快。其基本思想是通过用 h(x, θ) 替换 f(x, z, θ) 中的 z 来计算 Lu，忽略模型缺陷 ζ。因此，上层数据的对数似然函数为：
$$
L_u = \ln L_u = \sum_{i=1}^{N_u} p(y_nu | x_nu, \theta, \Sigma_y, \Sigma_z) = \sum_{i=1}^{N_u} \left(-\frac{1}{2} \ln(2\pi) - \frac{1}{2} \ln | \Sigma_y | - \frac{1}{2} [y_nu - f(x_nu, h(x_nu, \theta), \theta)]^T \Sigma_y^{-1} [y_nu - f(x_nu, h(x_nu, \theta)]\right)
$$
完整数据和下层数据的对数似然函数分别为：
$$
L_c = \ln L_c = \sum_{i=1}^{N_c} \left(\ln p(y_nc | x_nc, z_nc, \theta, \Sigma_y) + \ln p(z_nc | x_nc, \theta, \Sigma_z)\right) = \sum_{i=1}^{N_c} \left(-\frac{1}{2} \ln(2\pi) - \frac{1}{2} \ln | \Sigma_y | - \frac{1}{2} [y_nc - f(x_nc, z_nc, \theta)]^T \Sigma_y^{-1} [y_nc - f(x_nc, z_nc, \theta)]\right)
$$
$$
L_l = \ln L_l = \sum_{i=1}^{N_l} p(z_nl | x_nl, \theta, \Sigma_z) = \sum_{i=1}^{N_l} \left(-\frac{1}{2} \ln(2\pi) - \frac{1}{2} \ln | \Sigma_z | - \frac{1}{2} [z_nl - h(x_nl, \theta)]^T \Sigma_z^{-1} [z_nl - h(x_nl, \theta)]\right)
$$
然后可以通过最大化完整数据的对数似然来估计参数 θ、Σ_y 和 Σ_z：
$$
L_f = L_c + L_l + L_u
$$
例如，使用梯度下降算法。注意 Lu 关于 θ 的梯度为：
$$
\frac{\partial L_u}{\partial \theta} = -\sum_{i=1}^{N_u} \left[ y_nu - f(x_nu, z_ni, \theta) \right]^T \Sigma_y^{-1} \left(\frac{\partial f(x_nu, z_ni, \theta)}{\partial \theta} + \frac{\partial f(x_nu, z_ni, \theta)}{\partial z_n \cdot \partial h(x_nu, \theta)}\right) | z_n = h(x_nu, \theta)
$$
然而，需要注意的是，插值方法通常是误设的，因为它忽略了下层噪声 ζ。当下层函数 f(x, z, θ) 关于 z 是非线性的时，如果不使用完整数据或下层数据，这种误设会导致参数估计器的渐近偏差，如下命题所述。

命题 1
考虑在 (4) 和 (5) 中定义的混合模型。假设 f(x, z, θ) 对于某些 x 和 θ 至少是非线性的，并且 f 和 h 在 θ 和 z 上都是二次连续可微的。那么，如果只使用上层数据（即不使用完整数据或下层数据），插值估计器 θ 是渐近有偏的。证明见 B 部分。

2.3 不确定性量化
不确定性量化是建模的一个重要部分，因为它使我们能够估计预测的可靠性和变异性。在贝叶斯框架中，我们可以通过计算给定观测数据 D 的模型参数 θ 的后验分布来表征不确定性。这个后验可以围绕其模式 θ^ 进行近似（Gelman 等人，1995）：
$$
p(\theta | D, \Sigma_y, \Sigma_z) \approx N(\theta^, V)
$$
其中 θ^ 是在前一小节中得到的 θ 的估计值。协方差矩阵 V 由下式计算：
$$
V = [-\left(L_f\right)"(θ^)]^{-1}
$$
其中 L_f 是如方程 (14) 和 (17) 中定义的完整数据的对数似然，其海森矩阵由下式给出：
$$
(L_f)_{ij}"(θ) = \frac{\partial^2 L_f}{\partial \theta_i \partial \theta_j}
$$
这可以通过有限差分数值近似计算。注意，对于上层特定数据点，这个后验的计算需要在使用 MCEM 算法中进行重要性抽样。

给定一个新的输入 x*，我们用 z* 和 y* 分别表示下层输出和上层输出的相应预测。显然，z* 的预测均值是 z^*=h(x*,θ^)，y* 的预测均值是 y^*=f(x*,z^*,θ^)。现在我们考虑预测的不确定性。显然，参数的不确定性将传递到模型的预测不确定性。这可以通过给定 θ～N(θ^,V) 的正态近似来计算。通过对 θ^ 进行泰勒展开，我们有：
$$
h(x*,θ) \approx h(x*,θ^) + \nabla_\theta h(x*,θ^)(θ - θ^)
$$
因此，z* 的分布可以近似为均值为 z^*=h(x*,θ^) 和协方差矩阵为 Σ^z* 的正态分布。协方差矩阵是由参数协方差矩阵 V 和输出变量 Σ_z 的固有变异性组合得出的（Draper, 1995）：
$$
\S^z* = (\nabla_\theta h(x*,θ^))^T V(\nabla_\theta h(x*,θ^)) + \Sigma^z
$$
同样，y* 的分布可以近似为正态分布 (y^*,Σ^y*)。这里，协方差矩阵 Σ^y* 捕捉了参数和下层输出变量对上层输出变量的不确定性传播，计算如下：
$$
\S^y* = (\nabla_z,θf(x*,z^*,θ^))^T [Σ^模型的性能是通过均方根误差（RMSE）、平均绝对误差（MAE）和负对数预测密度（NLPD）指标来评估的，这些指标定义如下：

对于高层模型：
RMSE = 1/N × ∑_{i=1}^{N} (yi ? y^i)^2
MAE = 1/N × ∑_{i=1}^{N} |yi ? y^i|
NLPD = 1/N × ∑_{i=1}^{N} (?log(p(yi|y^i))

其中，yi 是模型计算出的真实值，y^i 是预测值，p(yi|y^i) 表示在给定预测值 y^i 的情况下观察值 yi 的预测密度。当预测服从均值为 y^i、方差为 σi^2 的正态分布时，p(yi|y^i) 对应于在 yi 处评估的正态概率密度函数的值。NLPD 考虑了预测区间（即不确定性）时的预测能力。如果预测值等于真实值且预测方差为零，NLPD 将达到最小值。低层模型的指标定义类似。

我们考虑了三种样本量（ny, nz）的情况：(49,23)、(36,36) 和 (16,64)，分别代表高层数据点多于、等于或少于低层数据的情景。

为了减轻随机变化的影响并提供可靠的模型性能估计，我们进行了10次独立实验，每次实验都是在给定的初始条件下随机生成样本点。RMSE、MAE、NLPD、估计参数（a, b, c）以及估计方差 Σ^y 和 Σ^z 的平均结果如表1所示。图1展示了每个实验情景中随机选择一个的预测结果。

表1. 基于十次重复实验的MCEM和Plug-in方法比较。

对于每种方法，表中还列出了参数 Σ^y、Σ^z、RMSE、MAE、NLPD、估计参数 a、b、c 以及估计的方差 Σ^y 和 Σ^z 的值。

我们考虑了三种样本量（ny, nz）的情况：(49,23)、(36,36) 和 (16,64)，分别代表高层数据点多于、等于或少于低层数据的情景。

为了减轻随机变化的影响并提供可靠的模型性能估计，我们进行了10次独立实验，每次实验都是在给定的初始条件下随机生成样本点。RMSE、MAE、NLPD、估计参数（a, b, c）以及估计方差 Σ^y 和 Σ^z 的平均结果如表1所示。图1展示了每个实验情景中随机选择一个的预测结果。

表1. 基于十次重复实验的MCEM和Plug-in方法比较。

在第二个模拟中，我们使用以下非线性模型生成数据：
(24) y = a1 + e^(-bx) + cz^2 + ?,
(25) z = b1 + x^2 + ?，
其中真实参数 a、b、c 分别为 [0.4, 1.2, 0.7]。误差项 ? 和 ? 是独立的，并且服从均值为零、方差分别为 Σ? = 0.0016 和 Σz = 0.0016 的正态分布。协变量值 {xyi} 和 {xzi} 都是在 [-6, 4] 范围内均匀分布的随机样本。

然后执行与模拟1中相同的模型识别和预测程序，结果总结在表2中。图2展示了每个实验情景中随机选择一个的预测结果。

表2. 基于十次重复实验的MCEM和Plug-in方法比较。

总体而言，两项模拟 study 的结果显示四种方法在性能上存在明显差异。在大多数情景中，MCEM 方法在参数估计准确性、预测准确性和不确定性量化方面取得了最佳性能。这体现在其通常较小的 RMSE 和 MAE 值上，更重要的是，其 NLPD 值始终较低，表明预测分布校准得更好。Plug-in 方法通常能提供有竞争力的点预测，并且在计算上更为简单，但其不确定性估计不太可靠，因为低层模型的缺陷没有完全传播到高层。因此，尽管在某些情况下其 RMSE 和 MAE 可能与 MCEM 接近，但其 NLPD 值通常更差。Two-step 方法通常比 Complete 方法表现更好，这是预期的，因为它至少利用了低层的具体数据来估计潜在的低层关系，然后再构建高层预测。然而，由于 Two-step 策略将低层和高层估计阶段分开，因此它没有充分利用混合结构，其预测和不确定性量化在大多数情况下仍不如 MCEM。

在四种方法中，Complete 方法的整体性能最弱。由于仅依赖完整数据并丢弃高层和低层的特定观测值，它未能充分利用可用信息，因此生成了不太准确的参数估计和较差的预测性能。

综上所述，这些结果突显了在混合框架内联合利用所有可用数据的优势。图1和图2还表明，在需要准确不确定性量化时，显式考虑潜在的低层不确定性是至关重要的。

在真实数据应用中，我们专注于 MCEM 和 Plug-in 方法，因为这两种方法是在所提出的统一框架内开发的，并且在模拟研究中都显示出比 Complete 和 Two-step 方法更好的性能。

在本节中，我们将 MCEM 和 Plug-in 方法应用于化学计量学中的一个真实世界数据集，具体是关于人类皮肤角质层（SC）的分配系数。众所周知，SC 是化学物质渗透皮肤的主要屏障。SC 由异质材料组成，主要包括脂质和角质细胞（由蛋白质（角蛋白）组成的死细胞）。化学物质在脂质和角蛋白之间的分配系数对于理解皮肤渗透性至关重要；它们是开发皮肤药代/毒代动力学机制模型所必需的。

本研究使用的数据集来自 Wang 等人（2010年）报告的实验测量结果。它包含了来自多个来源的异质观测数据，包括95对（Kow, Ksc）、30对（Kow, Kpro）和15对（Kow, Klip）。辛醇-水分配系数 Kow 作为主要输入变量，涵盖了广泛的化学疏水性范围。在模型识别之前，数据经过了仔细的一致性和单位一致性筛选。为了保持机制关系的物理可解释性，没有进行额外的标准化处理。

这个选定的示例为评估所提出的框架提供了一个有代表性的测试案例。尽管其规模适中，但该示例包括部分观测变量，并展示了在化学计量学应用中常见的多层结构，使其特别适合评估所提出的混合识别框架。

分配系数 Ksc 被建模为三个主要组分的函数：化学物质与蛋白质（角蛋白）、脂质和水之间的分配系数，如方程（26）所示：
(26) Ksc = φproKpro + φlipKlip + φw + ?sc。
此外，Kpro 和 Klip 分别通过以下公式建模：
(27) Kpro = ρproρwθ1Kowθ2 + ?pro,
(28) Klip = ρlipρwKowθ3 + ?lip
这里 Kow 是辛醇-水分配系数，可以很容易地获得大量化合物的数值（无论是实验上还是通过计算化学方法）。参数 θ1、θ2 和 θ3 是未知的，需要估计，ρpro、ρlip 和 ρw 分别代表蛋白质、脂质和水的体积密度，根据先前的研究分别为 ρpro = 1.37 g/cm3、ρlip = 0.9 g/cm3 和 ρw = 1 g/cm3（Nitsche, Wang, & Kasting, 2006）。误差（?sc、?pro、?lip）假设服从均值为零、方差分别为 Σ?sc、Σ?pro 和 Σ?lip 的正态分布。

需要注意的是，在我们的混合模型中，Ksc 对应于高层输出，而 Kpro 和 Klip 对应于低层输出。上述模型并不是独立的，因为 Ksc 取决于 Kpro 和 Klip。此外，每种模型可用的信息量也不同。对于这个特定的示例，从实验上分离 SC 中的蛋白质和脂质以测量 Kpro 和 Klip 是困难的。因此，如数据集中所述，报告的 Ksc 数据相对较多。目前，在这种情况下缺乏公认的模型识别方法。一种经验方法（Wang 等人，2010年）是直接从可用数据中识别模型（28），然后使用模型（26）中估计的 θ3，再使用 SC 分配数据来估计 θ1 和 θ2。显然，这种方法并不理想，因为它没有利用所有可用的数据来识别多个模型中出现的参数（例如，θ1 和 θ2 出现在模型（26）和（27）中，而 θ3 出现在模型（26）和（28）中）。显然，上述模型（26）、（27）和（28）在数学上很简单，因为它们本质上是单变量非线性回归问题，用于将 Kow 与性质 Ksc、Kpro 和 Klip 关联起来。可以认为，纯粹数据驱动的模型，如人工神经网络和高斯过程回归，可以直接从数据构建并取得良好的拟合（可能还有预测）结果。然而，在许多工程和应用化学领域，人们仍然更倾向于使用像这样的机制模型，因为它们提供了更清晰的物理解释。例如，模型（26）表明，可以从各个相（蛋白质、脂质和水）的分配中计算出 SC 中的分配。

为了利用所有可用数据，我们应用 MCEM 和 Plug-in 方法来识别模型参数，并通过留一法交叉验证（LOOCV）预测高层输出 Ksc 以及低层输出 Kpro 和 Klip。表3展示了平均 RMSE、MAE 和 NLPD 的结果，图3展示了预测的高层和低层输出以及观测数据。

表3. 通过留一法交叉验证得到的化学物质分配系数的平均 RMSE、MAE 和 NLPD。

从表3和图3可以看出，MCEM 方法在多个方面都优于 Plug-in 方法。具体来说，MCEM 方法的高层输出 Ksc 和低层输出 Kpro 和 Klip 的 RMSE 和 MAE 值都更小，表明预测准确性更好。此外，MCEM 方法的所有输出的 NLPD 值都显著低于 Plug-in 方法，表明 MCEM 提供了更准确和更精确的预测不确定性。图3还进一步证实了 Plug-in 方法高估了预测方差，正如在模拟示例中所看到的。

在本文中，我们介绍了一个新混合模型，该模型仅适用于数据集不完整的情况，并提供了一种统一的混合模型识别和不确定性量化方法。与需要所有变量完整观测值的传统方法不同，所提出的框架允许使用部分观测信息在多个模型层进行参数估计和不确定性传播。这种能力使得该方法在数据可用性有限的复杂建模问题中特别有用。我们提出了两种方法，即蒙特卡洛期望最大化（MCEM）算法和 Plug-in 方法，用于模型识别和不确定性量化，并在模拟研究中将它们与两种基准策略（Complete 方法和 Two-step 方法）进行了比较。此外，还使用涉及角质层中化学物质分配系数的真实数据应用来证明所提出框架的实际效用。此外，还建立了关于混合模型可识别性和 Plug-in 参数估计器渐近偏差的理论结果。这些结果提供了关于四种方法的优点和局限性的见解，以及它们在实际应用中的指导。从数值示例可以看出，四种方法表现出不同的优势和劣势。MCEM算法在预测不确定性量化方面始终表现出色，这从模拟和实际应用案例中较低的负对数预测密度（NLPD）值可以明显看出。这一优势归因于MCEM能够考虑潜在变量，并通过迭代采样来更全面地处理不确定性，从而改进参数估计。然而，这种优势的代价是计算复杂度增加，这在处理大型数据集或计算资源有限的情况下可能成为一个显著的缺点。值得注意的是，完全贝叶斯推断方法（如基于MCMC的方法）虽然提供了全面的不确定性表征，但通常伴随着较高的计算成本，特别是对于包含潜在变量和不完整数据的混合模型。相比之下，所提出的基于MCEM的框架在保持可靠不确定性量的同时， computation效率更高，这一点在模拟研究和实际数据应用中都得到了验证。这种权衡使得该方法在计算效率至关重要的实际化学计量建模场景中特别有吸引力。

另一方面，Plug-in方法在点预测准确性方面也能取得有竞争力的结果，并且由于其相对简单和高效的计算性能而具有吸引力。但对于较为简单的模型或对精确不确定性量化要求不高的场景，它的计算效率使其成为一个不错的选择。然而，Plug-in方法并不能完全将低层模型的不确定性传播到高层预测中，因此在模拟研究和实际数据示例中显示出较不可靠的预测不确定性量化结果。

两种基准策略——Complete方法和Two-step方法——进一步强调了充分利用不完整混合数据结构的重要性。Complete方法仅使用完整数据并忽略高层和低层的特定观测值，在四种方法中表现最差，这表明丢弃部分观测数据可能会导致信息丢失和估计及预测性能下降。Two-step方法首先估计底层模型，然后将预测的底层输出插入高层模型，通常比Complete方法表现得更好，因为它利用了低层的特定数据。但由于这两个阶段是分开处理的，Two-step方法仍然未能充分利用混合结构的联合特性，且不确定性传递效果不如MCEM。总体而言，这些比较突显了在所提出框架内联合建模所有可用数据的优势。

在实际应用中，选择哪种估计算法取决于可用数据集的结构和具体问题的需求。当有足够的观测数据时，可以直接进行参数估计。当某些变量未被观测到时，所提出的框架可以通过混合模型内的耦合关系进行间接识别。对于需要稳健不确定性量化的应用，尤其是在存在潜在变量或层次依赖性的情况下，MCEM表现更好。相反，如果主要目标是高效地进行点估计而不对预测不确定性有太高要求，Plug-in方法可以提供一个更简单且快速的替代方案。基准测试结果还表明，Complete方法和Two-step方法可以作为有用的参考策略，但它们通常不如所提出的方法有效，因为它们没有充分利用不完整的多层次数据结构。数值实验表明，在需要准确不确定性量化且计算资源允许的情况下，MCEM通常是更优的选择；而在计算资源受限的情况下，Plug-in方法更适合快速估计。

未来的研究可以探索结合Plug-in方法的计算效率和MCEM提供的稳健不确定性表征的混合方法，例如，先用Plug-in方法进行初始参数估计，然后通过有针对性的MCEM细化来平衡计算效率和预测准确性。在本研究中，MCEM中的模型参数采用非信息性先验，以使参数估计过程主要依赖于观测数据中的信息。虽然这种方法减少了主观先验假设可能引入的偏差，但当可用数据集有限时，可能会导致不确定性区间变宽。研究结果对不同先验规格（如信息量较少的先验）的敏感性可以提供进一步的见解，并为未来的研究指明方向。此外，将这些方法扩展到更复杂的模型（如贝叶斯神经网络或高斯过程）中，可以进一步提高它们在处理复杂实际问题时的适用性，其中不确定性在决策中起着关键作用。此外，所提出的框架主要适用于使用现有数据集进行离线混合模型识别。不过，该方法也可以扩展到在线或连续环境中，即随着时间的推移新观测数据不断出现的情况。在这种情况下，可以通过纳入新收集的数据并重新运行MCEM程序（更新似然评估）来迭代更新参数估计。开发高效的在线实现是未来研究的一个有趣方向。

**附加信息：**
作者声明没有利益冲突。

**资助：**
本研究部分由中国奖学金委员会（编号202008060052）和莱斯特大学（编号CSE-MA5-WANG）资助。

**CRediT作者贡献声明：**
Hongri Cong：方法论、形式分析、研究调查、数据整理、验证、初稿撰写。
Bo Wang：概念化、监督、方法论、审稿与编辑。
Tao Chen：概念化、监督、方法论、审稿与编辑。