做古生物形态序列的人大概都经历过这样一个阶段手头攒了十几条来自不同剖面、不同层位的时间序列每条序列的采样点不多单独拿一条去拟合随机游走模型参数的标准误大得没法看AIC 也拉不开差距。这时候最容易犯的错是把几条序列的均值一拼凑成一条长序列扔进去跑。这么干通常能得到一个看着很漂亮的趋势但这个趋势很可能只是不同剖面的测量基准、地层对比误差堆出来的假象。paleoTS 这个 R 包里的fitMult就是为这类场景准备的喂给它一个装着多条序列的列表再加一个模型名它把同一个简单模型逐条套上去把似然加总返回一个能横向比较的合并结果。下面这篇把我自己在实际项目里用 fitMult 的完整链路写下来——数据怎么进、模型怎么挑、结果怎么读、哪些地方会不声不响地出错。1. 先搞清楚 fitMult 到底合并了什么1.1 单序列拟合的天花板在哪里在聊 fitMult 之前得先承认一个现实古生物时间序列的数据质量天生就撑不起复杂的模型。一条序列能有 15 到 30 个采样点已经算丰富了而每个采样点上的均值本身还带着测量误差和个体变异。你拿这么一条序列去估一个三参数的随机游走模型自由度所剩无几参数的置信区间宽得能跑马。我在早期项目里做过一个统计二十来条形态序列单条拟合 GRWgeneral random walk时vstep的相对标准误普遍在 40% 以上个别序列甚至超过 100%这种结果拿去做速率比较基本没有意义。问题不在于模型错了而在于信息量不够。一条序列能提供的关于演化速率的证据被采样密度和样本量死死卡住。这时候人的本能反应是多凑点数据但古生物数据不像实验数据可以随时补采你手里有什么就得用什么。fitMult 的价值就在这儿它不去改变单条序列的信息量而是把多条序列的证据按似然相加的方式汇总起来让你至少能在这批序列整体上更符合哪种演化模式这个层面给出一个判断。1.2 合并的是似然不是参数这是全文最需要说清楚的一点也是我见过最多人误解的地方。fitMult 做的事情用一句话概括对列表里的每一条序列用同一个模型设定分别拟合然后把各自的 log-likelihood 加起来得到一个总的 logL、总的参数个数 K、总的样本量 n进而算出合并的 AIC、AICc、BIC。它并不会让多条序列共享一套参数。也就是说如果你跑了fitMult(y, model GRW)返回结果里vstep是一列数每条序列一个值而不是一个所有序列共用的值。这一点和分层贝叶斯模型里部分汇聚的思路完全不同。为什么 paleoTS 要这么设计因为在古生物场景下所有物种的演化速率相同这个假设通常站不住脚强行共享参数会把真实的生物学差异抹平但所有物种都遵循同一类演化模式比如都是随机游走而非停滞这个假设反而是有意义、值得检验的。fitMult 支持的是后者。理解了这一点你就知道它适合回答什么问题、不适合回答什么问题你想回答的问题fitMult 是否合适原因这批序列整体更像随机游走还是停滞合适合并似然后可直接比 AICc每条序列各自的演化速率是多少合适参数按序列分别给出所有序列的平均演化速率是多少需二次加工要自己对参数做加权汇总各序列的速率是否存在真实差异不合适它不做方差成分分解某个环境因子是否影响演化速率不合适需要回归型或分层模型提示拿到 fitMult 的返回对象后第一件事是用names()和str()把它拆开看一眼。不同版本的字段名可能有出入$parameters、$logL、$K、$AICc、$all.fits这几个是我常用的但不要凭记忆硬写字段名。1.3 什么情况下不该用它有三种情况我会明确避开 fitMult。第一种序列之间的时间跨度差了一个数量级。比如一条序列跨了 200 万年另一条只跨了 5 万年。这种情况下合并 logL 在数学上没问题但解释上很危险因为随机游走的vstep是每步的方差而步的定义依赖采样间隔跨度和采样密度差太多时合并出来的整体结论会被长序列主导。第二种序列的性状量纲不统一。一篇文章里同时放体长毫米和角度度各自的vstep单位完全不同合并 logL 只是在做数值加法但你在正文里没法把它们放在一张表里比较速率。真要做得先标准化而且标准化之后vstep的解释就变成了每步变化多少个标准差得在方法里写清楚。第三种你真正想问的是演化模式在类群之间是否不同。这时候正确做法是分组各跑一次 fitMult然后比较两组的合并 AICc而不是把所有序列一锅炖。2. 把剖面数据整理成 fitMult 认得的列表2.1 单个 paleoTS 对象的三个硬性要求fitMult 的第一个参数是一个 list里面每一项都得是合法的paleoTS对象。构造对象最常用的是as.paleoTS()核心参数是mm各采样点的性状均值、vv各点的方差通常是个体变异、nn各点的样本量、tt各点的时间或年龄。这里有三个坑几乎每个新手都会踩至少一个。时间必须从老到新也就是升序。如果你习惯用距今年代Ma BP那数值是从大到小排列的直接喂进去不会报错但所有基于相邻点差分的计算全反了。paleoTS 里对tt有检查但检查的严格程度取决于版本我见过不报错但结果荒谬的情况。稳妥做法是构造完对象后手动确认all(diff(x$tt) 0)必须返回 TRUE。单位统一。tt用 Ma 还是 kamm用毫米还是厘米全篇必须一致。我曾经在一个合作项目里把两条分别用厘米和毫米记录的序列合到一起跑vstep差了四个数量级图上完全看不出来因为画的是各自的相对变化。vv不能是 0。有些老数据集只给了均值没给方差如果vv全填 0似然函数里会出现除零或者方差为负的问题优化器会给出莫名其妙的收敛码。我的处理方式是给一个基于测量误差的下限比如 0.01 倍的均值平方并在方法里注明这是假设值。2.2 从一张长表拆成列表的写法实际数据通常是一张长表一列是序列 ID一列是时间一列是均值一列是方差一列是样本量。用split()按 ID 拆开再对每一组调用as.paleoTS()是最省事的路子library(paleoTS) # dat 是长表taxon / age / size / var / n dat - dat[order(dat$taxon, dat$age), ] # 先按 ID、再按年龄升序排 ts_list - lapply(split(dat, dat$taxon), function(d) { d - d[order(d$age), ] # 组内再排一次双保险 as.paleoTS(mm d$size, vv d$var, nn d$n, tt d$age) }) # 逐个自检 sapply(ts_list, function(x) all(diff(x$tt) 0)) sapply(ts_list, function(x) length(x$mm))注意split()出来的列表顺序是按 ID 的字母序不是你原始数据里的出现顺序。这个顺序在 fitMult 里会一路传递到结果的参数矩阵行名上所以后面整理结果表的时候最好用一个显式的 ID 向量对齐别指望行号。另外一个细节as.paleoTS()有个reset.time之类的参数不同版本命名不同作用是把时间平移到从 0 开始。对单条序列拟合来说平移只影响anc的取值不影响vstep和 logL但对于多条序列如果你不做平移各序列的anc就是各自时间轴上的绝对值横向比较没有意义。我一般会统一把所有序列的时间都平移到起点为 0这样anc就变成了起始状态比较好解释。2.3 长度不齐和样本量悬殊怎么办fitMult 允许列表里的序列长度不同这一点比很多人想的宽松。但它不会替你处理某条序列只有 4 个点这种情况。4 个点拟合一个三参数模型自由度只剩 1优化器给出的参数基本等于瞎猜而且这条序列的 logL 会以很奇怪的方式影响合并结果。我的经验阈值是单条序列至少 8 个采样点且至少一半的点上有nn 3的样本量。低于这个门槛的序列要么合并进相邻层位要么单独拿出来做描述性统计而不进 fitMult。样本量悬殊的问题更隐蔽——一条序列每个点 50 个个体另一条只有 2 个个体后者的vv估计极不稳定会通过似然函数放大噪声。处理办法是给vv设一个下限或者在nn很低的点上适当放大vv把不确定性显式地写进模型。注意不要为了让所有序列看起来整齐而做插值补齐。古生物序列的采样点对应的是真实地层位置插出来的点在生物学上不存在会人为压低vstep的估计。3. 模型菜单五种简单模型各自假设了什么3.1 参数个数与生物学含义paleoTS 里 fitMult 支持的模型名常用的就那几个字符串。选模型之前必须清楚每个模型在讲什么故事否则 AICc 选出来的东西你没法解释。模型名讲的故事主要参数参数个数 KGRW一般随机游走有方向性漂移anc, mstep, vstep3URW无偏随机游走无方向性漂移anc, vstep2Stasis停滞围绕某个均值波动theta, omega2StrictStasis严格停滞均值固定omega1Trend方向性趋势含时间趋势项anc, mstep, vstep3OU均值回归被拉向最优值anc, theta, omega, vstep4GRW和URW的区别只有一点前者允许每一步的期望变化不为零后者强制为零。这个区别在短序列上几乎测不出来所以很多论文里 GRW 和 URW 的 AICc 只差不到 2属于无法区分。Trend和GRW在数学形式上很接近区别在于趋势项的处理方式实际数据里两者结论常常一致我一般同时跑看哪个更简洁。Stasis是最容易被低估的模型。很多研究默认演化就是变化跑完 GRW 就收工但如果数据其实更符合停滞硬套 GRW 会把噪声解释成方向性漂移得出缓慢但持续演化的错误结论。我的习惯是每次分析必跑Stasis和StrictStasis作为基线哪怕最后没选它们。提示上表的 K 值是常见设定下的结果anc是否计入参数、StrictStasis的均值是否完全固定不同版本的处理略有差别。拿不准的时候直接看summary()输出里的 K或者用length(fit$parameters)核对。3.2 Joint 和 AD 两种似然写法的取舍fitMult 的method参数通常有Joint和AD两个选项这是很多人直接跳过、但其实影响很大的一个选择。Joint是联合似然它把整条序列当作一个多元正态分布来处理需要构造完整的协方差矩阵。这个矩阵的元素由随机游走或 OU的结构决定规模是采样点数的平方。好处是用上了全部信息包括不相邻点之间的关系坏处是对采样点多的序列计算量大而且协方差矩阵的估计对异常值敏感。AD是 ancestral-descendant 方法只看相邻两个采样点之间的差异把每一步当成独立增量来处理。它的计算快得多对长序列友好但丢掉了非相邻点之间的相关信息估计效率略低。我的实际选择标准很朴素序列平均长度小于 25 个点时用Joint超过就用AD如果两种方法给出的vstep差了一倍以上说明数据里有强自相关或者异常点得回头查数据而不是纠结方法。另外同一篇文章里所有模型必须用同一种 method混用会让 AICc 失去可比性。3.3 选模型的正确顺序我不建议一上来就把六个模型全跑一遍然后看 AICc 最小的那个。正确的顺序是先看数据的形态再定候选集最后比信息准则。第一步把每条序列画出来。看它是围绕一个水平线上下波动像 Stasis还是单调爬升像 Trend还是无规则游走像 GRW。这一步纯靠眼睛但能排除掉一半不可能的模型。第二步根据问题定候选集。如果研究对象是形态的长期稳定性候选集就是Stasis、StrictStasis、URW如果研究的是某个环境事件后的快速变化候选集是GRW、Trend、OU。第三步才是跑 fitMult 比 AICc。这时候比较的是 2 到 3 个模型而不是 6 个结论也更容易讲清楚。顺带说一句AICc 差值小于 2 的时候不要硬选直接说数据无法区分这几个模型这在审稿时反而更站得住脚。4. 跑通一次完整分析从列表到结果表4.1 最小可运行代码先用模拟数据把流程跑通确认环境没问题再上真实数据。paleoTS 自带几个模拟函数sim.GRW()、sim.Stasis()、sim.OU()这些用来造测试数据很方便library(paleoTS) set.seed(2024) # 造 6 条序列3 条随机游走3 条停滞 g1 - sim.GRW(ns 20, ms 0, vs 0.05) g2 - sim.GRW(ns 18, ms 0, vs 0.05) g3 - sim.GRW(ns 22, ms 0, vs 0.05) s1 - sim.Stasis(ns 20, omega 0.8) s2 - sim.Stasis(ns 19, omega 0.8) s3 - sim.Stasis(ns 21, omega 0.8) pool - list(g1, g2, g3, s1, s2, s3) names(pool) - c(g1, g2, g3, s1, s2, s3) # 同一个模型套到所有序列上 fit_grw - fitMult(pool, model GRW, method Joint) fit_sta - fitMult(pool, model Stasis, method Joint) summary(fit_grw) summary(fit_sta)这段代码里唯一需要留意的是names(pool)。给它命名之后结果的参数矩阵行名就是这些名字后面整理表格能省不少事。不命名也能跑但行名会变成默认的数字索引多条序列混在一起时极易对错行。4.2 结果对象怎么读summary()打印出来的信息量有限真要提取数值还得自己动手。我常用的几个字段names(fit_grw) # 先看有哪些字段 fit_grw$logL # 合并的 log-likelihood fit_grw$K # 合并的参数总数 fit_grw$AICc # 合并的 AICc fit_grw$parameters # 参数矩阵每行一条序列 fit_grw$all.fits # 各序列单独的拟合结果$parameters是最有价值的部分它把每条序列的anc、mstep、vstep排成矩阵。拿到之后第一件事是看vstep那一列的分布如果某条序列的vstep比其他序列大两个数量级多半是那条序列里有个别离群点或者采样间隔特别短。$all.fits是一个列表里面存着每条序列单独拟合的结果。这个字段常被忽略但它非常有用——你可以拿它跟$parameters对照检查合并拟合有没有把某条序列的参数带偏。理论上 fitMult 是逐条独立拟合的两者应该完全一致如果发现不一致说明模型设定或者数据有问题。合并的 AICc 怎么算出来的就是把各序列的 logL 加起来K 加起来n 加起来然后套用 AICc 的标准公式。这个做法成立的前提是各序列相互独立。来自同一剖面的相邻序列可能不独立共享地层对比误差这种情况在解释时要打折扣或者考虑把这几条序列先合并。4.3 把结果整理成可用的表格跑完三四个模型之后手动抄数字太容易出错。写一个小函数统一提取tidy_fit - function(fit, label) { data.frame( model label, logL as.numeric(fit$logL), K as.numeric(fit$K), AICc as.numeric(fit$AICc), n_seq if (!is.null(fit$parameters)) nrow(fit$parameters) else NA ) } res - rbind( tidy_fit(fit_grw, GRW), tidy_fit(fit_sta, Stasis) ) res$dAICc - res$AICc - min(res$AICc) res[order(res$dAICc), ]这段代码里as.numeric()是保险起见。有些版本返回的是带名字的向量直接塞进 data.frame 会变成多列看着莫名其妙。n_seq那一行的判断也是同样道理防止字段缺失时整个函数报错。出来的表加上dAICc一列按差值排序基本就是论文里模型比较表的样子。记得在表注里写清楚 K 是合并值、n 是总采样点数、用的是哪种 method这三条信息不写审稿人一定会问。5. 我踩过的坑和对应的排查链路5.1 时间顺序反了不报错但全错这是最阴的一个坑。数据表里如果时间列是距今年代数值从大到小而你直接按原始顺序构造对象tt就是递减的。某些版本的as.paleoTS()会给出警告某些版本什么都不说。跑出来的vstep看着正常logL也正常但mstep的符号会整体反过来趋势的方向完全颠倒。排查链路是这样的第一步diff(x$tt)看是不是全正第二步如果发现有负值不要只对那一条排序检查是不是整个数据表的排序逻辑有问题第三步重新构造对象后把新旧两次拟合的mstep对比一下符号应该翻转。我现在的习惯是在构造列表的循环里直接断言stopifnot(all(diff(d$age) 0))一行代码省掉后面几天的返工。5.2 收敛失败与离谱的参数fitMult 内部的优化用的是数值优化器遇到病态数据会不收敛。表现形式通常是vstep跑到 1e-8 或者 1e6 这种极端值或者mstep的标准误是 NA。这时候不要急着换模型先诊断数据。我的排查顺序是先看这条序列的采样间隔分布如果相邻点间隔差异极大比如有的点隔 1 万年有的隔 50 万年随机游走的每步假设就不成立得考虑把间隔作为权重或者干脆重采样到近似等间隔。再看vv有没有异常小的值vv接近 0 会让似然函数在边界上打转。最后看序列长度少于 6 个点的序列基本别指望能稳定估出三个参数。如果数据本身没问题还是不收敛可以试着手动给初值。paleoTS 的拟合函数一般支持传初值虽然 fitMult 的接口里不一定直接暴露但你可以先用fitSimple()对单条序列试跑找到稳定的初值再回到 fitMult。5.3 AICc 不可比的几种隐蔽情况同一批序列上跑不同模型AICc 是可以比的。但下面这几种情况会让比较失效而且不会有任何报错提示。第一种method不一致。Joint和AD算出来的 logL 定义不同混着比等于拿苹果比橘子。第二种序列集合不一致。比如跑 GRW 的时候列表里有 12 条序列跑 Stasis 的时候因为某条序列报错被剔掉了变成 11 条。这时候 AICc 的差值毫无意义。我的做法是跑之前先固定列表任何一条出问题就整体停下修数据不要跳过继续。第三种vv的处理不一致。有的模型在拟合过程中会重新估计vv有的会固定住这会导致 logL 的口径不同。稳妥做法是在所有模型里统一固定vv为观测到的个体方差并在方法里写明。5.4 标准误不等于显著性fitMult 的结果里通常带参数的标准误。很多人在表里给vstep标上星号说显著大于零。这个推断在古生物时间序列里要非常小心因为vstep的抽样分布严重偏斜而且参数在边界附近接近 0时标准的 Wald 型置信区间根本不可靠。我自己的做法是不用星号直接报参数和标准误如果要检验速率是否显著大于零用似然比检验比较vstep自由估计和固定为接近零的模型或者用参数自举。这两种方法在计算上都不复杂但结论可靠得多。6. 把它用顺手的几个进阶做法6.1 一次跑完所有候选模型手动一个个调 fitMult 太笨。写个循环把候选模型一次性跑完顺便把结果收成一张表run_models - function(pool, models c(GRW, URW, Stasis, StrictStasis, Trend)) { out - lapply(models, function(m) { fit - tryCatch( fitMult(pool, model m, method Joint), error function(e) NULL ) if (is.null(fit)) return(NULL) data.frame(model m, logL as.numeric(fit$logL), K as.numeric(fit$K), AICc as.numeric(fit$AICc)) }) do.call(rbind, out) } tab - run_models(pool) tab$dAICc - tab$AICc - min(tab$AICc) tab[order(tab$dAICc), ]tryCatch()这一步很关键。真实数据里总有某条序列在某个模型下不收敛没有这层保护整个循环会中断前面跑的全白费。返回 NULL 之后表格里少一行你能一眼看出来哪个模型有问题。6.2 单条序列的结论和总体结论打架怎么办这是实际分析里最常见的困境合并结果显示 GRW 明显优于 Stasis但翻看$all.fits会发现有三四条序列单独看其实更像停滞。这时候不要简单地以总体为准。我的处理方式是把序列按单条拟合的 AICc 偏好分成两组然后对两组分别跑 fitMult比较两组的合并结论。如果两组确实呈现不同的模式那这本身就是一个值得报告的发现——说明这批序列里存在模式异质性而不是统一遵循某一种演化方式。这种先分后合的思路比强行给一个总体结论要有说服力得多。反过来如果分组之后每组内部还是很乱那就说明单条序列的信息量根本不足以支撑模型判别这时候老老实实报告数据无法区分这几种模式是最诚实的选择。6.3 把拟合结果叠到序列图上纯看表格很难形成直觉把拟合曲线画到原始序列上会清楚很多。单条序列的拟合对象一般有对应的plot()方法可以直接叠par(mfrow c(2, 3)) for (i in seq_along(pool)) { plot(pool[[i]], main names(pool)[i]) # 若 all.fits 中有对应单条拟合可在此叠加拟合线 }六条序列放一张图上哪条像游走、哪条像停滞一眼就能看出来。这个图我一般放在补充材料里正文只放合并结果的表两者配合读者既能看整体结论也能自己判断每条序列的情况。最后分享一个我在实际操作中的体会fitMult 这种同模型套多条序列的工具最大的价值不在于给出一个漂亮的 AICc 排名而在于逼你把每条序列都检查一遍。我做过的好几个项目里真正改变结论的往往不是模型选择而是在整理列表的过程中发现某条序列的地层对比有问题、或者某个层位的样本量记录有误。跑代码只要几分钟查数据要几天但这几天的价值远大于那几分钟。