《Methods in Ecology and Evolution》：Early warning signal for river-borne diseases with almost no data
1 引言
水传播疾病，特别是河流性疾病，对环境和医疗保健造成巨大负担。每年约有340万人死于与水相关的疾病，病原体包括大肠杆菌、霍乱（cholera）和沙门氏菌（salmonella）。病原体的分布因人类运输和全球变暖而迅速变化。早期预测病原体存在有助于减缓疾病传播，但这需要分析从微观宿主-病原体水平到宿主间相互作用的多种变量，以及生态、社会、经济和人口因素。然而，在疫情爆发初期，由于监测和诊断测试的限制，实验数据往往稀疏，而环境协变量数据通常详细。此外，与许多向所有方向均匀扩散的入侵物种不同，河流中的病原体通过特定的空间结构传播，该结构可用有向无环图（DAG）表示。
这些挑战限制了预测模型的使用。机制模型如偏微分方程能够解释复杂的空间过程，但难以拟合到异质性结构的河流网络，且预测准确性低。标准的易感感染移除模型（SIR）及其空间扩展在捕捉疾病传播方面有成功记录，但强调结构方面而非协变量的影响，可能未充分利用详细的环境数据。简单统计模型如广义线性模型（GLM）无法捕捉变量间的复杂交互，而高级机器学习模型如长短期记忆（LSTM）需要比可用数据更多的数据。河流网络文献引入了线性模型的成功扩展，称为混合模型，通过纳入空间自相关随机效应来捕捉河流结构和协变量效应。然而，该模型在疫情早期需要大量数据实例。
隐马尔可夫模型（HMM）用于模拟未观测变量（隐藏状态）的演化，该变量随时间“发射”可观测的发射（如环境观测）。作为图形概率模型，HMM能够模拟具有未知动力学的过程，并处理训练和测试时的缺失值。然而，它们无法模拟比路径更复杂的空间过程。HMM已被扩展到具有每个节点多个出边的树结构、多个入边和DAG、有序树和平行链。但这些适用于河流的扩展HMM假设，当多条河流合并时，下游河流的污染可能性等于上游合并河流污染可能性的平均值，这在生物学上不合理，因为如果任何支流的相邻上游区域被污染，下游区域很可能也被污染。此外，尽管河流主要由汇入的支流组成，但它们也可能分支形成分流和岛屿，这些未被树结构捕捉。更重要的是，基于HMM的模型都不允许发射之间的互连，即它们假设环境观测在给定隐藏状态条件下是独立的，这可能显著降低模型的预测性能。
本文扩展了HMM以处理这些问题。我们允许隐藏状态通过任意DAG连接，父节点通过Noisy-OR对子节点执行影响。此外，通过使用树增强朴素贝叶斯（TAN），我们允许每个切片内的发射直接相互依赖。作为案例研究，我们模拟了加拿大阿尔伯塔省奥尔德曼河鳟鱼旋转病的传播。我们检查了模型预测准确性随可用测试结果数量从0增加到100的变化，并将结果与混合模型进行比较。
鳟鱼旋转病由引起，这是一种鲑科鱼类（包括鳟鱼、鲑鱼和白鱼）的微观寄生虫。幼鱼易受旋转病影响，死亡率高达90%。这种疾病对鱼类种群和渔业构成威胁，由于对休闲钓鱼和水产养殖的不利影响，可能造成经济损失。1950年代首次在美国记录，旋转病在美国西部和东北部各州流行。2016年8月，加拿大首次在班夫国家公园约翰逊湖和阿尔伯塔省鲍河流域的鱼类中确认旋转病的存在。2017年4月，加拿大食品检验局宣布鲍河和奥尔德曼河流域感染旋转病，阿尔伯塔省其余地区被指定为缓冲区。
2 方法
2.1 模型
考虑任意河流并使用网格将其划分为方形像素。河流可以通过DAG建模，其中节点是像素，边表示水流。具体地，如果水从像素i流向像素j且两个方形像素通过角或边相邻，则像素i链接到像素j。相应地，如果存在从该节点到节点j的有向路径，则该节点是节点j的上游像素，下游像素相反。像素的河流是其所有上游和下游像素的集合。对于每个像素i，定义响应变量（隐藏状态）Yi，指示特定病原体的存在（1）或缺失（0）。我们通常将病原体存在称为污染。还考虑一些告知病原体存在的环境观测，例如河流温度，堆叠在实值向量Xi中，表示像素i处的值。目标是为每个像素找到给定该像素河流所有环境观测的病原体存在概率。
我们通过假设变量如何相互依赖来接近这个问题。让Pa(Yi)表示节点i的直接上游像素的隐藏状态变量。第一个假设是，给定隐藏变量Yi的父节点，它变得独立于其余上游响应变量。这本质上是一个马尔可夫假设：每个像素的污染可能性直接仅取决于其相邻上游像素，而不是更上游的河流。第二个假设是，给定响应变量Yi，环境观测Xi变得独立于所有其他环境观测和隐藏变量。即，如果我们知道像素i的污染状态，那么知道其他像素的污染状态和环境观测不会改善我们对像素i环境观测的知识。
上述两个假设导致所有变量的联合概率分解为转移概率分布和条件观测概率分布（也称为发射概率分布）的乘积。转移概率分布确定病原体从直接上游像素传输到像素本身的概率，条件观测概率分布给定像素污染状态的环境观测似然。
我们如下模拟转移概率。假设病原体从任何污染父节点（直接上游像素）独立地以概率p传输到下游像素。那么下游像素清洁（未污染）的概率是病原体未从任何污染父节点传输的概率，即(1-p)m，其中m是节点i的污染直接上游像素数。如果所有父节点清洁，我们仍分配污染概率p0给外部因素。这产生Noisy-OR模型。
关于观测概率分布，传统方法是做出“朴素”假设，即观测在给定响应变量条件下独立，导致朴素贝叶斯（NB）结构。我们通过允许每个观测除了隐藏状态外最多依赖于一个其他观测来放松这一假设，产生TAN模型。TAN可以可视化为隐藏状态和环境观测上的图：对于每个Xi(k)，从条件变量绘制入边。结果图是TAN结构。TAN参数是概率项。我们模拟P(Xi| Yi)而不是P(Yi| Xi)以简化参数。使用贝叶斯规则和上述方程，可以回答任何概率查询。
如果每个节点恰好有一个父节点，且观测概率采用NB结构，模型将成为众所周知的HMM。然而，响应变量在HMM中不可观测，称为隐藏状态，但发射观测，因此称为发射。我们的模型因此将HMM扩展到具有TAN结构观测的部分可观测DAG-Markov模型，其中发射被与节点疾病状态相关的环境观测取代。环境观测在概念上不同于发射的经典概念，因为它们不是由疾病引起而是概率相关。
尽管TAN-HMM未直接将河流流量或其他环境观测纳入转移概率，但通过发射概率间接这样做。这是因为水文信息被纳入有助于发射概率的环境观测中。具体地，像素i给定观测的污染概率通过对联合概率分布中所有其他响应变量边缘化获得。该方程包括发射概率项P(Xi| Yi)，如果作为环境观测包含，则编码水文信息。
给定具有已知参数的扩展HMM，我们推广前向-后向算法以获得每个像素给定其河流环境观测的污染概率。然而，必须首先使用期望最大化（EM）等算法从数据学习参数。数据通常仅限于部分环境观测且无或少数隐藏状态变量。环境观测中的缺失值可以在训练模型时使用EM与隐藏变量一起估计或插补。学习称为监督（如果响应变量在某些像素已知）或无监督（否则）。我们假设河流网络已知且未学习它，但从数据学习TAN结构。
2.2 设置和数据
奥尔德曼河的鳟鱼旋转病通过取M. cerebralis寄生虫的存在作为响应变量建模。我们使用阿尔伯塔环境和公园收集的数据集，因此本研究不需要伦理批准。代码在Zenodo可用，使用代码的分步手册在支持信息中提供。我们使用200 m × 200 m像素大小的网格研究奥尔德曼河流域。我们仅考虑与河流重叠的像素并分配网格代码。表示河流结构的DAG有根节点和叶节点。每个节点最多有四个父节点和五个子节点。部分节点只有一个子节点，部分节点只有一个父节点。从2017年和2018年实验室进行的鱼类和蠕虫测试中，响应变量Yi的值可用于113个像素，仅占总像素数的0.07%。其中，30个实例报告阳性，其余83个阴性。对应每个像素i，我们考虑九个环境观测，产生观测向量Xi。数据集的每个实例由网格代码标记，如果不缺失，包含九个观测和响应变量的值。九个观测中的四个包括缺失值，三个（电导率、水pH和溶解氧）的值在98%的数据实例中缺失。所有九个观测是连续变量。
2.3 学习和评估
我们比较了四个模型在预测未知测试结果像素响应变量方面的性能：混合模型、无河流结构的TAN模型、具有河流结构和TAN观测的HMM（扩展TAN-HMM；主要模型）和具有河流结构和NB观测的HMM（扩展NB-HMM）。我们考虑了10批实验：
批次实例1。每批从随机选择113个数据实例中的2个用于模型训练开始，其余111个用于测试。我们确保两个中一个是阳性测试结果，另一个是阴性，因为混合模型需要两者执行分类。
批次实例2。对于每批的下一个实例，从111个中随机选择3个更多实例，导致总共5个训练和108个测试实例。
批次实例3。然后随机添加5个更多，导致总共10个训练和103个测试实例。
批次实例4-10。从那里开始，每次迭代添加10个新实例，直到获得总共100个训练和13个测试实例。
零训练实例的情况也执行，但仅一次，因为它不涉及任何随机化。
所有环境变量首先离散化为三个等宽箱（‘低’、‘中’、‘高’），然后使用bnstruct包中的k-最近邻（k-NN）插补一次。k-NN搜索仅在环境变量空间执行；未使用空间信息（例如坐标或河流网络距离）。产生的离散化和插补数据集保持固定，并用于在所有批次中拟合和评估所有四个模型以确保公平比较。
TAN的结构是“离线”学习的，即一次从隐藏状态和环境变量的数据样本中学习，独立于河流DAG，使用R中的bnlearn包。这产生隐藏状态Yi和发射Xi上的相同固定树增强图。然后将结构纳入网络以学习TAN的参数以及转移概率和根像素（无父节点）的污染概率。参数在每批中分别使用该批的训练数据集学习。每批的训练数据集包括所有164,985个实例的环境数据以及该批中包含的污染状态（响应变量）。参数通过最大化模型给定训练数据集的似然估计。没有已知的易处理算法找到全局最大值；然而，可以使用Baum-Welch算法高效导出局部最大似然，这是著名EM算法的特例。算法迭代地找到参数的最大似然估计，分两步。在期望步中，使用前向-后向算法计算每个像素响应变量污染期望概率，给定初始参数和每个像素的环境变量。该信息然后在最大化步中用于更新参数，最大化观测数据似然。此过程继续直到参数估计收敛。对于每批，我们初始化转移参数为p = 0.5和p0= 0.5，根节点响应变量为P(Yi) = 0.5。我们设置初始发射概率分布为均匀分布。
对于无网络河流的TAN，我们使用bnlearn包从训练数据集学习参数。对于学习具有NB结构发射的扩展HMM，我们开发了TreeHMM包，公开可用。
混合模型使用方程预测响应变量，其中响应变量向量、环境观测矩阵、空间自相关随机变量向量、参数向量和误差向量。向量捕获网络河流中像素之间的欧几里得距离自相关、流连接像素之间的自相关（尾上）以及流连接和流未连接像素之间的自相关（尾下）。为了检查混合模型，我们首先使用R包raster、sf和openSTARS将旋转病测试结果和环境观测数据转换为空间河流网络（SSN）对象。然后使用SSN包训练混合模型。
最后，我们使用接收者操作特征曲线下面积（AUC）分数在测试数据集上测试每个模型的性能，平均超过所有10批。
3 结果
扩展TAN-HMM在仅两个像素的疾病测试结果可用时得分0.7 AUC。混合模型在此情况下表现0.57 AUC，当测试结果增加到5时约0.55 AUC。对于训练大小2-10，扩展TAN-HMM优于混合模型；在几乎所有训练大小上也优于扩展NB-HMM和仅TAN模型。从20个测试结果开始，混合模型优于扩展TAN-HMM，并在提供100个测试结果时达到接近完美的性能。扩展TAN-HMM在此点表现0.9 AUC。
对于无疾病测试结果可用的特定情况，扩展TAN-HMM基于环境变量将像素分类为两类。如果任一类正确匹配病原体存在或缺失，分类结果将为0.71 AUC。混合模型无测试结果时不进行预测。
我们检查了不同模型组件如何贡献于预测准确性。如果忽略河流连接，导致响应变量无空间相关的模型，平均AUC将降至0.6（当有2-5个测试结果可用时），但降至0.5（更多测试结果时），与随机猜测一样差。如果忽略环境变量之间的依赖性，导致简单NB结构，性能下降但不严重。
学习的TAN结构被检查以回答一般查询，形式为找到任意像素给定该像素单个环境变量的污染概率。病原体存在可能性随水pH、干扰、电导率和气温增加，但随坡度、水质和流量减少。
4 讨论
模拟河流性疾病具有挑战性，由于非均匀结构、通常少数可用测试结果和未知动力学。如果期望早期预警信号，挑战加剧，通常接近无测试结果可用。我们推广HMM以处理奥尔德曼河旋转病数据的这些问题。我们将河流划分为像素，每个像素引起病原体的存在由二元响应变量模拟，并从每个像素链接到其相邻较低海拔像素，产生DAG。我们通过TAN模拟每个像素的环境观测。尽管旋转病可用测试结果少且一些环境变量在爆发前多年测量，学习的扩展TAN-HMM能够以0.7 AUC预测河流上疾病状态，比最先进混合模型高约0.13 AUC。混合模型在20个测试结果可用后开始优于。结果因此表明扩展HMM可以在少数到几乎零测试结果可用时作为早期预警信号，一旦疾病通过初始出现状态并收集足够测试结果，可能被混合模型取代。
预测性能在通过DAG互连响应变量后显著增加，强调河流结构和邻近效应在病原体传播中的作用。实际上，与无任何空间连接的TAN相比，其表现与随机分类器一样差，允许每个像素响应变量依赖于其上河流响应变量将AUC增加0.12（少数可用测试结果）和约0.4（所有100个可用时）。注意尽管连接响应变量的DAG是单向的，信息流是双向的。例如，像素的污染可能性可以使用其下游像素之一估计，反之亦然。在响应变量的马尔可夫假设可能被违反的应用中——例如，当病原体传播依赖于宿主扩散时——宿主存在可视为隐藏状态，而病原体存在连同其他特征作为发射。然后EM算法仍可用于估计发射值包括病原体存在。
有趣的是，即使扩展HMM无旋转病测试结果学习，仅使用环境观测的空间模式，仍能提供研究区域的二元分类，匹配疾病测试结果0.7 AUC。这是在假设两个二元值正确映射到病原体存在和缺失的情况下；否则，AUC翻转到0.3 AUC。尽管有些令人惊讶，结果与机器学习中无监督学习文献一致。这表明河流中一些区域是除M. cerebralis外其他病原体的潜在栖息地，受相同环境变量影响。此外，我们使用的大多数环境变量，如气温和坡度，应可供监测管理人员访问。因此，为了检查先前清洁河流中其他病原体的存在，管理人员可以通过向扩展HMM提供相关环境观测来优化调查工作，以识别最可能污染的区域。
重要的是澄清，虽然我们已确定模型为早期预警系统，其主要功能可能更准确地描述为一般识别高风险区域，而非特别精确定位有旋转病的网格单元。这种区别突出模型作为一般风险预警工具而非早期检测机制的潜力，特别是给定其早期分类性能与爆发数据量无关的一致性。然而，即使少数早期侵染点的存在特别有助于验证高风险区域已正确识别。在此背景下，模型最有力地应用于入侵早期阶段。
更好纳入空间相关性的能力可能解释混合模型在训练大小20个实例或更多时优于扩展HMM。这与偏差-方差权衡一致：扩展TAN-HMM将空间依赖性限制到每个像素的直接上游节点（父节点），并将环境变量之间的依赖性限制到树结构。因此，该模型有偏差，因为它无法捕捉一些结构；然而，它具有较低的模型复杂性，这意味着其参数的近最优值可以用有限数据估计，即低方差。相比之下，混合模型利用空间自相关随机效应，允许它在更多数据可用时更稳健地捕捉河流不同部分之间的固有空间依赖性。因此，混合模型可能有较低偏差但较高方差，因为它需要更多数据来适当估计其参数。
给定单个变量条件，模型表明当水质低或电导率高时，污染概率超过0.5。这些可能性是边际的，即当未提供其他环境测量或污染数据时适用。这些和其他数据驱动模式与关于疾病如何与坡度、流量、水质和人为干扰相关的文献一致。此外，将气温作为河流温度的近似，较高温度似乎为奥尔德曼河寄生虫生长提供更合适的基质，这在麦迪逊河也有报告。关于水特性，结果与先前研究相似，表明较高电导率增加污染可能性。水pH在单变量图中可能显得特别信息丰富；然而，由于它仅对1%–2%像素观测且其余从其他协变量插补，此类显示仅诊断性且应谨慎解释。
扩展TAN-HMM的AUC在0和20个训练标签之间几乎平坦（0.70）。我们相信这是因为无监督组件——环境协变量耦合河流网络——已经良好排名风险；前几个标签主要设置类别极性并略微校准参数，因此基于排名的判别（AUC）变化很小，直到模型看到更多标签。此外，尽管早期AUC变化适中，地图重新校准：随着训练大小增加，高概率区域的空间范围和强度转移。
生态上，早期AUC平台表明环境过滤和流动连通性在初期入侵期间已经构建风险表面；少数早期检测主要验证和定向该表面而非重塑它。这与基于生态位推理平行但与环境生态位模型对比，后者从先验范围推断物种环境属性——通常使用GARP等算法——然后将这些属性投影到新景观以预测潜在建立。相比之下，我们的模型从当前范围协变量和河流网络原位学习类似生态位的风险表面，无需先前范围的任何信息。因此，我们的模型最适合没有先前范围长期研究历史的新引入。
我们的模型为生态学中其他应用开辟了新途径。它有潜力应用于引起河流性疾病的其他病原体，包括沙门氏菌、霍乱（霍乱弧菌）和鮰爱德华菌，以及入侵物种，如斑马贻贝。通常，此类疾病的数据在环境变量方面详细，但在测试结果方面极其有限，这是我们的模型表现最佳的情况。此外，污染测试结果可以由0到1之间的数字表示以捕捉混合污染状态，例如当同一像素有两个不同测试结果可用时。这是因为EM学习算法利用充分统计量，意味着它分配概率给二元响应变量中的缺失值，而不是将它们插补为确定的0或1。该模型还可能预测水生数量变化和水灾，如洪水，每个像素的洪水水平作为响应变量，地理、水文和气候特征作为环境变量。地震预测是另一个应用，因为地震断层遵循类似河流的结构，可以划分为像素，地震发生作为响应变量，地震特征作为环境变量。更广泛的应用包括癌症通过血管传播，其中血管被像素化，响应变量是循环肿瘤DNA的存在，环境变量可以是微血管密度和细胞游离DNA或液体活检。
我们离线学习TAN结构，即在纳入河流网络之前；然后我们学习TAN参数。可以扩展模型以也基于河流网络学习TAN结构。此外，可以考虑使用其他模型，如线性回归或神经网络，而不是TAN或NB。这允许环境变量之间更复杂的关系。它还允许我们检查贝叶斯网络以外的模型是否受益于结构化响应变量。另一个限制是使用离散数据。实现连续分布变量的算法，如高斯分布，留作未来工作。
作者贡献
Pouria Ramazi、Russell Greiner和Mark A. Lewis概念化并监督研究，并贡献于分析。Prajwal Bende和Arezoo Haratian实现代码并执行分析。Pouria Ramazi领导写作，所有作者贡献于写作。
利益冲突声明
作者无利益冲突。
数据可用性声明
数据可通过
https://doi.org/10.5281/zenodo.17468107获取。
