paleoTS fitMult：多时间序列合并拟合与AICc模型选择-育师

1. 为什么要把多个时间序列捆在一起拟合

手上攒了四五条形态测量序列，分别来自不同剖面、不同层位，每条序列的采样点数量、时间跨度、个体数都不一样。你想回答的问题其实很朴素：这些序列背后是不是同一种演化模式在起作用？如果一条一条地用fitSimple去拟合，你会得到四五个互不相干的 AICc，接下来要么凭感觉挑一个，要么把结果硬凑成一张表，审稿人一句“这些 AICc 怎么加起来的”就能把你问住。paleoTS包里的fitMult就是为这种场景准备的：它把多个paleoTS对象放进一个列表，对每条序列拟合相同结构的简单模型（GRW、URW、Stasis 等），再把对数似然相加，给出一个合并的 AICc。你拿到的不再是一堆散点，而是一个可以横向比较的总体指标。

这个函数在古生物形态演化、微体化石尺寸序列、生物地层学里的性状变化研究里用得非常频繁。它的门槛不高，但坑不少：pool参数一改，模型含义就完全变了；序列独立性假设一旦被打破，似然相加就是自欺欺人；AICc 的分母还会在你序列很短的时候突然变成负数。这篇东西就是把这些细节一次讲透，从数据整理到结果报告，给你一套能直接抄作业的流程。

1.1 fitMult到底解决了什么痛点

单序列拟合的核心问题是信息量太薄。一条十几个采样点的序列，拟合一个三参数的 GRW 模型，残差自由度所剩无几，AICc 的惩罚项会大得离谱，模型排名基本靠运气。fitMult的思路是把多条序列的似然拼起来：每条序列仍然各自估计自己的参数，但模型结构必须一致，最后把 logLik 加总。这样总样本量上去了，AICc 的惩罚相对更合理，模型之间的区分度也就出来了。

这里要说清楚一个容易误解的点：fitMult默认不是让所有序列共享一组参数。它做的是“同结构、异参数”——每条序列有自己的祖先均值、自己的步长方差，只是大家都用“随机游走”这个骨架。真正共享参数要靠pool = TRUE，那是另一回事，后面单独讲。很多人第一次用fitMult，以为它自动做了参数池化，结果解释结果时把“各序列参数不同”说成了“共同参数”，这是典型的概念混淆。

从实现上看，fitMult内部就是对你传入的列表做一次循环，对每个元素调用单序列的拟合函数，再把logLik、K、n累加。所以它的结果对象里保留了$allFits，你可以随时把每条序列的单独拟合结果拎出来看。这个设计很实用：合并指标用来做模型选择，单序列结果用来画图、看参数、做后续分析。

1.2 和单序列fitSimple的本质区别

fitSimple接收一个paleoTS对象，返回一个paleoTSfit；fitMult接收一个列表，返回的也是paleoTSfit，但里面的logLik、AICc、K、n都是合并后的。字段名一样，含义完全不同，这是新手最容易栽跟头的地方。你如果拿fitMult的结果去调plot或者提取$par，会发现结构对不上——因为$par在pool = FALSE时并不是一个扁平向量，而是藏在$allFits里，每条序列一份。

另一个区别是模型选择的粒度。fitSimple回答的是“这一条序列更像随机游走还是更像停滞”；fitMult回答的是“这一组序列整体上更像随机游走还是更像停滞”。前者是个案判断，后者是组间判断。写文章的时候，如果你的论点是“这个类群在多个剖面上都表现为定向演化”，那fitMult的合并证据比四五条单序列各自显著要有力得多。

还有一点值得注意：fitSimple支持fitModels一次性拟合多个模型并直接比较，fitMult没有对应的批量接口。你得自己对每个模型调一次fitMult，再用compareModels把结果拼起来。多写几行代码而已，但流程上要心里有数。

1.3 什么时候不该用fitMult

fitMult的数学基础是序列间独立，似然才能直接相加。如果你的几条序列来自同一条剖面、共享同一批个体，或者时间上有重叠、地层上连续，那它们就不独立，合并似然会高估信息量，AICc 会偏向复杂模型。这种情形下要么把序列合并成一条长序列（如果时间轴能接上），要么在模型里显式考虑序列间的相关性，而不是硬套fitMult。

第二种不该用的场景是量纲或测量尺度不一致。比如一条序列测的是臼齿长度（毫米），另一条测的是体长（厘米），你还开pool = TRUE，那共享的vstep就是个没有物理意义的混合体。即便pool = FALSE，合并 AICc 也只在“模型结构可比”的意义上成立，解释时要小心。

第三种是序列太短。每条序列至少要有五六个采样点，fitMult才有意义。两三个点的序列拟合随机游走，参数比数据还多，结果全是噪声。与其硬凑，不如把短序列剔除或者合并。

2. 动手前的数据准备与对象构造

2.1 安装与版本确认

paleoTS在 CRAN 上有稳定版本，直接装就行。我习惯在脚本开头把版本号打出来，因为不同版本之间fitMult的参数名和返回字段偶有调整，尤其是pool的默认值和hess的支持情况。你复现别人结果翻车，十有八九是版本对不上。

install.packages("paleoTS") library(paleoTS) packageVersion("paleoTS")

如果你要用开发版的新特性，可以从代码托管平台装，但生产脚本里我不建议这么干，版本漂移会让你几个月后自己都跑不通。装完之后用?fitMult看一眼当前环境的函数签名，比对着任何教程都靠谱。

2.2 把野外数据整理成paleoTS对象

真实数据一般是表格，每行一个采样点，列包括序列编号、时间、均值、方差、个体数。整理成列表的通用模板如下：

library(paleoTS) df <- read.csv("measurements.csv", stringsAsFactors = FALSE) df <- df[order(df$series_id, df$time), ] split_list <- split(df, df$series_id) ts_list <- lapply(split_list, function(d) { as.paleoTS( mm = d$mean, vv = d$variance, tt = d$time, ni = d$n ) })

几个硬性要求：tt必须按升序排好，paleoTS不会替你排序；vv不能有零或负值，因为内部要做对数变换；ni至少为 1，缺了就用 1 填，但要知道这会改变抽样误差的估计。如果你的表里只有均值和标准差，记得先把标准差平方再填进vv。

构造完先别急着拟合，用summary(ts_list[[1]])看一眼对象长什么样。paleoTS对象的核心字段是$mm（均值）、$vv（方差）、$tt（时间）、$ni（个体数），还有一些派生字段。确认这几个字段的长度一致、没有 NA，再往下走。

2.3 序列质量检查清单

我在正式拟合前一定会跑一遍下面这张检查表，能挡掉后面八成的报错。

对应代码：

check_ts <- function(x) { c( n = length(x$mm), tmin = min(x$tt), tmax = max(x$tt), vv_min = min(x$vv), sorted = all(diff(x$tt) > 0) ) } sapply(ts_list, check_ts)

这张表跑完，你对手上数据的“健康状况”就有数了。哪条序列时间跨度特别短、哪条方差特别大，后面解释 AICc 的时候都用得上。

3. fitMult的核心参数与调用方式

3.1 函数签名逐项拆解

fitMult的签名在不同版本里略有差异，但核心参数稳定：

fitMult( ts.list, model = c("GRW", "URW", "RW", "Stasis", "Stasis+strict"), pool = FALSE, method = "BFGS", hess = FALSE, ... )

逐项说明：

model的可选值对应paleoTS的简单模型族。GRW是广义随机游走，含祖先均值、步长方差、每步漂移三个参数；URW是无偏随机游走，只有祖先均值和步长方差；Stasis是停滞模型，含最优值和约束强度两个参数；Stasis+strict是严格停滞，只有一个最优值参数。RW在多数版本里和URW参数个数相同，具体以你环境里fitSimple返回的$K为准，不要凭记忆写死。

3.2 pool参数：独立参数还是共享参数

这是fitMult最容易被误读的参数，值得单独拎出来讲。

pool = FALSE（默认）：每条序列独立估计自己的一套参数，模型结构相同，似然相加。总参数个数K = S × K_single，S是序列数。这相当于说“大家都用随机游走这个骨架，但每条序列的步长和起点可以不同”。适合序列来自不同地点、不同物种、时间尺度不完全可比的情形。

pool = TRUE：所有序列共享同一组参数，总参数个数K = K_single。这相当于说“这些序列背后是同一个演化过程，起点、步长、漂移都该一样”。这个假设非常强，只有在序列确实来自同一总体、时间尺度可直接对比时才成立。它带来的好处是参数少、AICc 惩罚轻，但代价是一旦假设不成立，模型会严重失真。

我的习惯是两种都跑，看 AICc 差多少。如果pool = TRUE的 AICc 明显更低，说明数据支持“共同过程”，这时候再报告共享参数才有底气；如果pool = FALSE明显更好，说明序列间异质性大，老老实实报各序列独立参数。两种结果都报出来，审稿人反而觉得你想得周全。

3.3 模型选择与参数个数K的核对

不同模型的参数个数不一样，pool一开，总数还会变。下面这张表帮你快速核对：

| 模型 | 单序列参数 |pool = FALSE总 K |pool = TRUE总 K | | GRW | anc, vstep, drift | 3S | 3 | | URW | anc, vstep | 2S | 2 | | RW | anc, vstep（以实际为准） | 2S | 2 | | Stasis | theta, omega | 2S | 2 | | Stasis+strict | theta | S | 1 |

S为序列数。拟合完直接用fit$K核对，别自己心算。我见过有人在文章里写“GRW 模型共 3 个参数”，但他明明拟合了 6 条序列、开了pool = FALSE，实际 K 是 18，审稿人一算就对不上。这种低级错误一旦被抓到，整篇文章的可信度都要打折。

4. 完整实操：从三个序列到一张AICc表

4.1 构造示例数据

没有现成数据的时候，用模拟数据把流程跑通最省事。paleoTS自带sim.GRW、sim.Stasis等函数，参数名如果和你本地版本对不上，用?sim.GRW查一下。

set.seed(2024) ts1 <- sim.GRW(ns = 15, ms = 0.02, vs = 0.05, vp = 0.02) ts2 <- sim.GRW(ns = 12, ms = 0.02, vs = 0.05, vp = 0.02) ts3 <- sim.GRW(ns = 18, ms = 0.02, vs = 0.05, vp = 0.02) ts_list <- list(A = ts1, B = ts2, C = ts3)

三条序列的采样点数量不同（15、12、18），这正好模拟了真实数据里长度不齐的情况。命名列表不是必须的，但强烈建议，因为后面提取$allFits的时候，名字能帮你对上号。

4.2 单序列试跑

正式合并之前，先挑一条序列用fitSimple试一下，确认数据能被正常拟合：

fit_single <- fitSimple(ts1, model = "URW") fit_single

看三个东西：$logLik是不是有限的负数，$convergence是不是 0，$K是不是符合预期。如果logLik是-Inf或者convergence非零，先别往下走，回去查数据。单序列都跑不通，合并只会把问题放大。

4.3 多模型批量拟合

确认单序列没问题，就可以对列表逐个模型调用fitMult：

fit_urw <- fitMult(ts_list, model = "URW") fit_grw <- fitMult(ts_list, model = "GRW") fit_stasis <- fitMult(ts_list, model = "Stasis") compareModels(fit_urw, fit_grw, fit_stasis)

输出会给出每个模型的 logLik、K、AICc、ΔAICc 和 Akaike 权重。compareModels接受多个paleoTSfit对象，fitMult返回的正好是这个类，所以能直接拼表。

如果要试pool = TRUE，再加一行：

fit_grw_pool <- fitMult(ts_list, model = "GRW", pool = TRUE)

然后单独比较fit_grw和fit_grw_pool，看共享参数假设是否被数据支持。

4.4 结果解读与模型平均

拿到 AICc 表之后，别急着只看第一名。我的解读顺序是：

先看 ΔAICc。第一名和第二名差小于 2，说明数据区分不开，两个模型都该报；差在 4 到 7 之间，算中等支持；大于 10，才可以放心说“明显更好”。再看 Akaike 权重，权重接近的时候做模型平均比硬选一个更稳妥。

模型平均没有现成函数，手动按权重加权参数就行：

w <- compareModels(fit_urw, fit_grw)$Akaike.wt # 假设第一行是 URW，第二行是 GRW vstep_avg <- w[1] * fit_urw$par["vstep"] + w[2] * fit_grw$par["vstep"]

注意pool = FALSE时$par不是扁平向量，参数在$allFits里，加权要按序列分别做。这也是为什么我建议先看清对象结构再写代码。

4.5 可视化与结果输出

原始序列的点图用plot(ts1)就能出来，横轴时间、纵轴均值、误差棒来自vv和ni。合并拟合结果本身没有现成的plot方法，我的做法是把每条序列的单拟合结果取出来，叠加到原始图上：

str(fit_urw$allFits, max.level = 1) pars_A <- fit_urw$allFits$A$par

导出结果表给合作者或者放进补充材料：

res <- data.frame( model = c("URW", "GRW", "Stasis"), logLik = c(fit_urw$logLik, fit_grw$logLik, fit_stasis$logLik), K = c(fit_urw$K, fit_grw$K, fit_stasis$K), AICc = c(fit_urw$AICc, fit_grw$AICc, fit_stasis$AICc) ) res$dAICc <- res$AICc - min(res$AICc) write.csv(res, "model_comparison.csv", row.names = FALSE)

这张表就是你文章里 Table 2 的雏形。列名建议用logLik、K、AICc、Delta_AICc、Akaike_weight，别自创缩写，读者一看就懂。

5. 常见问题与排查实录

5.1 收敛失败与初值问题

最典型的症状是convergence不等于 0，或者logLik返回NA。原因通常有三类：数据本身有异常值或零方差，初值离真实值太远，优化器选得不合适。

先查数据，用第 2 节的检查表过一遍。数据没问题的话，换优化器：

fit_try <- fitMult(ts_list, model = "GRW", method = "Nelder-Mead")

还不行就对数据做缩放。把每条序列的mm除以该序列的标准差，拟合完再把参数换算回去。缩放不改变模型结构，但能显著改善数值稳定性，尤其是序列量纲差异大的时候。

5.2 负数方差与参数越界

paleoTS内部对vstep一类参数做了对数变换，正常情况不会出现负估计。如果你在$allFits里看到负的vstep，八成是优化过程跑飞了，或者数据里有极端离群点。先画图看有没有离谱的采样点，有就剔除或降权，别硬拟合。

另一种情况是drift估计得特别大，导致预测均值飞出数据范围。这通常是序列太短、时间跨度太窄，模型无法约束漂移项。对策是退回URW，把drift固定为 0，用更简单的模型说话。

5.3 AICc为NA或K大于n

AICc 的公式里有n - K - 1做分母，当总样本量n不比参数个数K大多少时，分母趋近于零甚至为负，AICc 就会变成NA或负值。这在序列短、模型复杂、又开了pool = FALSE的时候特别常见。

判断方法很简单：fit$n和fit$K一比就知道。n - K - 1 <= 0就必须降复杂度。要么把GRW换成URW，要么把几条短序列合并，要么开pool = TRUE减少参数。没有别的捷径。

5.4 序列长度差异巨大时的权重问题

似然相加意味着长序列对总似然的贡献更大。一条 50 个点的序列和一条 8 个点的序列放在一起，合并结果基本由前者主导，后者相当于搭便车。这不是 bug，是似然加总的固有性质，但解读时要心里有数。

我的处理办法是：拟合前先看各序列的n，差距超过五倍的，要么单独报告短序列的结果，要么在讨论里明确说明“合并结果主要由长序列驱动”。如果短序列在生物学上很重要，可以考虑给它加权，但fitMult不直接支持权重，得自己改似然函数，一般不折腾。

5.5 常见报错速查表

6. 我踩过的坑与实操心得

6.1 不要盲目相信AICc排名

AICc 是相对指标，只在候选模型集合内部有意义。你如果只拟合了URW和Stasis，结论是“URW 更好”，但这不代表 URW 就是对的，只能说你没试过的模型里可能有更合适的。我的习惯是至少放URW、GRW、Stasis三个进去，GRW负责捕捉定向趋势，URW负责纯随机游走，Stasis负责停滞，三个骨架覆盖了最常见的演化模式。

另外一个经验是：序列数量少于 4 条时，AICc 的排名很不稳定，换个随机种子模拟数据，排名可能就翻过来。这种时候与其纠结第一名，不如把 ΔAICc 和权重都报出来，让读者自己判断证据强度。

6.2 参数提取与二次分析

pool = FALSE时，参数藏在$allFits里，每条序列一个paleoTSfit。提取步长方差的代码：

vstep_each <- sapply(fit_urw$allFits, function(f) f$par["vstep"])

pool = TRUE时，参数直接在$par里，一条序列一份。两种结构的差异一定要在脚本里注释清楚，否则几个月后你自己都会搞混。我一般会在拟合后立刻写一个小的提取函数，把结果整理成数据框，后面画图、做统计都方便。

还有一个细节：$allFits里的单序列拟合结果，logLik加起来应该等于合并的$logLik（在数值误差范围内）。这个恒等式可以用来验证结果有没有被意外修改。我习惯在关键步骤后跑一次核对，几秒钟的事，能省掉很多排查时间。

6.3 报告结果时的规范写法

文章里写方法部分，至少要说清楚四件事：用了哪些序列、每条序列多少采样点、拟合了哪些模型、pool开没开。结果部分给一张表，列logLik、K、AICc、Delta_AICc、Akaike_weight，正文里点明“ΔAICc 小于 2 的模型被认为具有相当支持”。

如果pool = TRUE和pool = FALSE的结果都报了，要明确说哪个是主分析、哪个是敏感性检验。我通常把pool = FALSE作为主分析，因为它假设更弱；pool = TRUE作为补充，用来检验“共同过程”是否站得住脚。这样写，审稿人基本挑不出毛病。

最后分享一个我自己用了很久的小习惯：每次跑完fitMult，把sessionInfo()和packageVersion("paleoTS")一起存进结果文件夹。古生物时间序列分析的可复现性，一半靠数据整理，一半靠环境记录。等到半年后合作者问你“当时那个 AICc 怎么算出来的”，你能直接把环境甩过去，比任何解释都管用。

paleoTS fitMult：多时间序列合并拟合与AICc模型选择