基于累积量的近似方法，实现物种分布的快速高效预测

《Ecological Modelling》：Cumulant-based approximation for fast and efficient prediction for species distribution

　　小森修|西草雄介|江口信人|久保田康弘
日本东京都武藏野市吉祥寺北町3-3-1，诚庆大学科学技术系，邮编180-8633

摘要
物种分布模型被广泛用于根据环境变量估算栖息地适宜性。在各种方法中，Maxent和泊松点过程是高效且常用的方法，但在处理大规模背景数据或高分辨率全球分析时，其计算成本会变得极高。当需依次对多种物种进行分析时，这一问题更为严重。为克服这一限制，我们提出了一种基于累积量近似方法的计算效率更高的物种分布模型，该方法应用于γ散度的损失函数。同时，我们还引入了一种带有L1惩罚的序贯估计算法，用于选择与物种分布密切相关的环境因素。基于累积量近似方法推导出的正则化几何平均法既具有较高的计算效率，又能保证估计精度。此外，通过将累积量近似方法应用于Maxent，我们证明了在环境变量线性组合服从正态分布的假设下，Maxent与费希尔线性判别分析是等价的。这一理论等价性使得估计过程可以极为快速地进行。通过对来自国家生态分析与综合中心的226种物种以及709种日本维管植物的两个数据集进行的模拟研究和分析表明，所提出的方法大幅降低了计算时间，同时保持了与Maxent相当的预测性能。R语言中的CBA包提供了用于模拟和实证分析的一些实现代码，便于在大规模物种分布模型中应用。重要的是，这种等价性还意味着可以使用glmnet包和高斯模型直接实现Maxent的二次近似，从而实现简单且高度可扩展的计算。

引言
物种分布模型被广泛用于研究物种与其相关环境特征之间的关系，从而阐明物种分布随时间变化的原因和机制（Elith和Leathwick，2009）。作为生态学、气候变化和保护生物学中的重要工具，当建模者与决策者之间能够良好协作时，这些模型可为保护决策提供支持（Guisan等人，2005；Guisan等人，2013）。用于构建物种分布模型的主要数据是仅记录物种存在情况的数据，这类数据缺乏物种不存在的信息，且可从博物馆和植物标本馆数据库中以相对较低的成本获取。相比之下，同时记录物种存在与不存在情况的数据则难以获得，因为其需要耗费大量时间和精力进行系统调查。在针对仅记录物种存在情况数据的物种分布模型中，Maxent（Phillips等人，2006）是一种信息论方法，可被用于多种用途，如保护规划以及研究入侵物种和气候变化的影响（Elith等人，2011）。Maxent旨在基于最大熵原理逐步最大化似然值（Dudík等人，2004），根据国家生态分析与综合中心包含的来自全球六个地区的226种物种的广泛数据，Maxent的性能位居前列（Elith等人，2006）。近期研究证实，与随机森林、XGBoost、支持向量机等其他机器学习方法相比，Maxent和提升回归树仍具有最佳性能（Valavi等人，2022）。理论研究表明，Maxent的估计结果与泊松点过程模型以及逻辑回归模型是等价的（Renner和Warton，2013；Warton和Shepherd，2010）。此外，Fithian和Hastie（2013）证明，无限加权逻辑回归也与泊松点过程模型等价。有关更详细的综述，可参考Komori等人（2023）的研究。近年来，物种分布模型已被应用于大规模数据集，以估算全球物种分布（Heshmati等人，2019；Lee等人，2021）。当使用高分辨率背景数据时，计算成本会变得极高（Phillips和Dudík，2008），通常通过随机抽样背景数据来缓解这一问题。Maxent中随机抽样的默认样本量为10,000（Hijmans等人，2024a），但有些研究认为，要获得可靠的结果，需要86,000个数据点（Renner等人，2015），而最新研究则建议根据ROC曲线下面积的值，使用50,000个数据点（Valavi等人，2022）。实际应用中，可通过不断增加背景数据点的数量，直到模型拟合和预测性能稳定为止（Phillips和Dudík，2008）。

本文提出了一种新的方法，通过利用所有背景数据点，并将累积量近似方法应用于Maxent的归一化常数或泊松点过程模型的强度函数数值积分，从而提高Maxent以及与之等价的泊松点过程模型的计算效率。这种基于玻尔兹曼机器学习中Plefka展开的Thouless–Anderson–Palmer近似方法（Kappen和de Borja Rodríguez Ortiz，1997；Tanaka，1998），还揭示了Maxent与费希尔线性判别分析之间的有趣关系：在正态分布假设下，Maxent参数的估计系数与费希尔线性判别分析的系数相等，这使得物种分布的估计速度极快。这一正态分布假设与对数高斯Cox过程相一致（M?ller等人，1998）。此外，我们将累积量近似方法应用于γ散度的损失函数（Basu等人，1998；Fujisawa和Eguchi，2008），由此衍生出一种名为正则化几何平均法的新方法，该方法兼具低计算成本和高预测精度。rGM的损失函数包含基于背景数据样本方差的正则化项，这有助于避免过拟合，优于基于几何平均散度得到的原始几何平均法。

本文的结构如下：首先，我们介绍泊松点过程模型及其γ损失函数，阐明它们与Maxent的关系及其统计特性。接着，我们将累积量近似方法应用于γ损失函数，证明在正态分布假设下，Maxent（即泊松点过程模型）与费希尔线性判别分析是等价的。随后，我们说明γ损失函数的一阶累积量近似对应于几何平均法，二阶累积量近似对应于正则化几何平均法。同时，我们还制定了rGM的估计算法，并在其中加入了用于变量选择的L1惩罚项。最后，我们在多种模拟场景和真实数据分析中评估了所提方法的性能，包括国家生态分析与综合中心的数据集（Valavi等人，2022）以及日本维管植物的数据集（Kubota等人，2015）。文章最后探讨了所提方法与现有方法的关系，并提出了进一步研究的建议。我们还从渐近阶数角度讨论了计算效率问题。CBA包中包含了用于模拟研究和真实数据分析的R语言代码，补充文件中也提供了如何使用CBA的简要说明。此外，附录中还给出了示例代码，展示了如何使用glmnet包高效实现Maxent的二次近似。

片段内容
泊松点过程及其快速估计算法
在泊松点过程回归模型中，我们处理的是仅记录物种存在情况的数据，这类数据由m个相互独立的地点{x1,…,xm}组成，同时还存在背景数据（也称为伪缺失数据），即研究区域A内的{xm+1,…,xn}这些地点。对于目标物种，我们考虑一个由斜率向量α和截距c参数化的强度函数λ0(xi,α,c)。通常，我们会采用对数线性模型，即logλ0(xi,α,c)=c+α?f(xi)，其中f(xi)是一个p维特征向量。

模拟研究
我们基于Cox过程设计模拟场景，其中强度λ(x,α)是随机变化的。设f(x)=(f1(x),…,fp(x))?表示一个p维特征向量。我们假设该特征向量的各分量之间的依赖结构由一个p×p相关矩阵R=(rjj′)定义，其元素满足rjj′=ρ|j?j′|(j,j′=1,…,p)。接下来我们考虑三种情况：
1. 高斯情况：我们生成f(xi)～N(0,R)；
2. 均匀情况：我们使用高斯Copula函数生成依赖变量fj(xi)～U(0,1)。

模拟研究
图1展示了在ρ=0.5、p=10、m=500、n=10,000的参数设置下，高斯情况、均匀情况和泊松情况下的模拟结果。在高斯情况下，根据推论1和附录E中的讨论，Maxent、rGM和费希尔方法的估计精度是有保障的。几何平均法的性能较差，因为它无法考虑特征向量之间的协方差。在计算效率方面，费希尔方法的计算速度约为Maxent的1000倍，比几何平均法和正则化几何平均法快约100倍。

讨论
在贝叶斯推断中，变分贝叶斯方法通过最大化变分自由能来近似后验分布，从而提高计算效率，这一方法基于变分密度（Friston等人，2007）。在均值场方法中，由于假设参数相互独立，该密度会被简化，同时还会运用拉普拉斯近似。后者为后验分布提供了局部高斯近似，有助于计算边际似然值。

CRediT作者贡献声明
小森修：撰写——初稿、软件、方法论、资金获取、概念构思。西草雄介：方法论、研究实施。江口信人：方法论、概念构思。久保田康弘：项目管理、资金获取。

生成式AI使用声明
作者们仅将生成式AI工具（如OpenAI的ChatGPT）用于文本编辑和提升写作清晰度。科学结果、讨论内容或结论中均未包含任何由AI生成的内容、想法或分析结果。作者对本文的内容负全责。

利益冲突声明
作者声明存在以下可能被视为潜在利益冲突的财务利益/个人关系：久保田康弘表示获得了日本文部省的财政支持；小森修表示获得了日本学术振兴会的财政支持。如有其他作者，他们也声明自己没有已知的可能影响本研究结果的财务利益或个人关系。

致谢
本研究得到了日本文部科学省“海洋生物大数据利用技术发展”项目（资助编号JPMXD152146988）以及日本学术振兴会的KAKENHI项目（资助编号22K11938）的支持。