SaTScan时空聚集检测:扫描统计、蒙特卡洛与七模型实践
简介:这份《SaTScan 中文说明书》面向公共卫生、流行病学、空间统计及地理信息分析的学习者与从业者,旨在降低英文文档门槛,帮助读者理解并操作 SaTScan 开展空间、时间及时空扫描统计,完成疾病聚集探测、随机性检验和异常信号评估。文档以单个 doc 文件交付,压缩包约 485KB,篇幅精炼,便于检索和随查随用。内容围绕软件用途、数据类型与统计模型展开,涉及离散与连续扫描数据、泊松模型、伯努利模型、时空重排模型、多项式与序模型、指数模型和正态模型,并说明地理坐标、背景人口、协变量、时间趋势及缺失数据的处理思路,同时给出软件下载安装、示例数据集和参数文件的使用提示。已有 2816 人学习,适合需要在疾病监测、犯罪学、生态学等空间分析场景中快速建立方法框架、理解输出结果的读者参考。
1. SaTScan 要回答的不是「有没有病」,而是「这堆病例扎堆扎得随机吗」
SaTScan 这类工具最容易被误解成“预测哪里会爆发疾病”。它真正在做的事要克制得多:给定一批已经发生的病例,以及它们的位置、时间,判断这些点在空间、时间或时空上是否偏离了随机分布。关键词是“偏离”和“随机分布”——必须先把随机分布定义成一个可计算的零模型,才谈得上偏离。这份中文说明书把这件事拆成两个层面:离散扫描统计假设位置是固定的、随机的是落在每个位置上的病例数;连续扫描统计假设位置本身就是随机的、随机的是研究区域内的点过程。搞懂这个区分,后面所有模型的取舍就都有了解释的锚点。
2. 似然比与蒙特卡洛:一个圆柱窗口怎么被判定为「显著聚集」
扫描统计量的形式化定义其实就一句话:在所有候选窗口里,找那个让“窗口内和窗口外各有一个独立的率”相对于“全局只有一个率”的似然比最大的窗口。这句话拆开有三个必须落地的工程问题——候选窗口怎么枚举、似然比怎么算、最大值分布怎么给 p 值。这三件事是 SaTScan 这类工具真正吃算力的地方,也是理解后续所有模型差异的基础。
2.1 候选窗口怎么枚举:圆形、圆柱与最大窗口比例
空间维度上,窗口是一个圆;时空维度上,窗口是一个圆柱——底面是地面上的圆,高是时间轴上一段连续区间。圆柱这个形状有明确动机:圆形让“邻近”有一个与方向无关的定义,而空间区间和时间区间的可分离性让窗口内的期望病例数可以写成空间部分与时间部分的乘积,计算量从二维降到一维。
枚举方式是暴力但可控的:圆心固定在一组候选位置上(病例位置或用户指定的格点),半径从小到大逐级放大;时间上起点和长度也逐级枚举。两个约束决定了枚举的边界:最大空间窗口比例和最大时间窗口比例,它们限制的是窗口最多覆盖多少比例的风险人口或研究期。常见取法是空间不超过总人口的 50%、时间不超过总研究期的 50%。这个限制是必要的,否则最大似然比总是倾向于把整个研究区装进窗口,得到的“聚集”就是全局本身。
还有一个认知偏差是:圆形窗口在空间上各向同性,如果真实聚集是一条沿河谷或公路的带状区域,单个圆覆盖不了,需要枚举多个重叠的圆去拼。椭圆窗口的引入就是为了贴合这种长条形聚集——两轴长度不同,扫描时分别枚举。
2.2 对数似然比:窗口内外两组计数在比较什么
先把离散泊松模型说清楚,因为其他模型基本是同一套路换掉分布假设。设研究区总病例数为 C,位置 i 的风险人口为 P_i,总人口为 P。零假设下病例按人口比例分布,位置 i 的期望病例数 μ_i = C · P_i / P。对任意候选窗口 Z,记窗口内观测病例数为 c、期望病例数为 μ,对数似然比为:
import math def poisson_llr(c, mu, C): """ c : 窗口内观测病例数 mu : 窗口内期望病例数(按背景人口比例分摊) C : 研究区总病例数 只检测"偏高"的聚集:c <= mu 时 LLR 置 0 """ if mu <= 0 or c <= mu or c >= C: return 0.0 inside = c * math.log(c / mu) # 窗口内:用自己的率替代全局率 outside = (C - c) * math.log((C - c) / (C - mu)) # 窗口外:剩余部分付出的代价 return inside + outside def scan_statistic(counts, expected, C): """在所有候选窗口上取最大 LLR,返回窗口下标与统计量值""" best_llr, best_idx = 0.0, -1 for i, (c, mu) in enumerate(zip(counts, expected)): llr = poisson_llr(c, mu, C) if llr > best_llr: best_llr, best_idx = llr, i return best_idx, best_llr逐项看这个表达式:第一项是窗口内部用“窗口自己的率 c/μ”替代“全局率 1”带来的似然增益;第二项是窗口外剩余部分在同样替代下付出的代价。因为全局总量守恒,两项此消彼长,只有窗口内的率确实高出全局、且高出的量足够大,总 LLR 才为正。这就是“随便挑一个病例数少的位置就说聚集”会被自动惩罚的原因——代价项不答应。
参数上有两个容易写错的边界:c >= C 或 mu <= 0 要直接返回 0,否则对数里会出现零或负数;c 和 mu 必须是同一时间窗、同一空间窗内的一致口径,混用会导致 LLR 虚高。实际源码里还会加一个数值下限,避免 log 在极端比例下溢出。
2.3 蒙特卡洛重排:p 值为什么是数出来的不是算出来的
扫描统计量是“所有窗口里最大的那个似然比”,这个最大值的解析分布对任意形状的窗口并不可得。SaTScan 用的是蒙特卡洛:先固定住数据的总量(比如总病例数 C),在零假设下重新生成一批数据,对每一批重算最大 LLR,然后看有多少次重算结果不小于真实观测值。
import random def monte_carlo_pvalue(obs_llr, pop_weights, C, windows, B=999, seed=7): """ obs_llr : 真实数据上的最大对数似然比 pop_weights: 各位置风险人口,用来按比例分配病例 C : 总病例数 windows : 每个候选窗口包含的位置下标列表,以及对应的期望病例数 mu B : 置换次数,决定 p 值的分辨率,最小可到 1/(B+1) """ rng = random.Random(seed) n_loc = len(pop_weights) exceed = 0 for _ in range(B): # 零假设:总病例数固定,按人口比例多项式重抽样到各位置 sim_counts = [0] * n_loc for loc in rng.choices(range(n_loc), weights=pop_weights, k=C): sim_counts[loc] += 1 # 用同一批候选窗口重算最大 LLR sim_max = 0.0 for locs, mu in windows: c = sum(sim_counts[i] for i in locs) sim_max = max(sim_max, poisson_llr(c, mu, C)) if sim_max >= obs_llr: exceed += 1 return (exceed + 1) / (B + 1)这段代码里最重要的不是循环,而是“固定总量、按比例重分配”这个约束:它保证了零假设下每个位置的期望病例数不变,变的是病例在位置间的随机波动。加 1 的分子分母是为了避免 p 值为 0,也是置换检验的标准写法。
这里有个必须记住的结论:p 值分辨率由置换次数决定。B=999 时最小 p 约 0.001,要报更小的 p 就得加次数,耗时线性增长。时空重排模型稍微特殊:窗口内的期望值是确定性算出来的(边际乘积分),置换的是时间标签,不是重新抽样病例数,但 p 值的数法不变。
3. 七个模型的取舍:同一个似然比框架,喂进去的数据长得不一样
如果把扫描统计量看成一个骨架,模型就是给这个骨架喂不同分布假设的“血肉”。SaTScan 提供泊松、伯努利、时空重排、多项、有序、指数、正态七种模型,看上去很多,实际是在回答同一个问题:你手里有什么数据,零假设下什么是随机的。
3.1 泊松与伯努利:有背景人口和没有背景人口,差别在分母
泊松模型的随机变量是每个位置上的病例计数,风险人口是外部给定的分母。零假设下病例数与该位置人口成比例。它适合数据能汇总成“某区域若干病例 + 一个背景人口数”的场景,比如按县、按邮区汇总的发病率数据。
伯努利模型的随机变量是个体的 0/1 状态:每个个体要么是病例,要么是对照。窗口内的分母是该窗口内个体总数,而不是一个抽象的风险人口。它适合病例-对照研究,比如出生缺陷研究里病例是缺陷新生儿、对照是同期无缺陷新生儿。两者的根本差别在于分母是否独立于病例计数:泊松模型里人口是输入、不随病例数变化;伯努利模型里病例和对照是同一批个体的两种状态,分母和病例数绑在一起。
这个差别直接体现在输入文件上:泊松模型需要人口文件(.pop),伯努利模型需要控制文件(.ctl)来记录每个位置有多少对照。混用了这两种模型,会出现“总病例数对得上但期望值全偏”的现象,程序不会报错,只是结果不可信。
3.2 时空重排模型:只拿到病例时,用边缘乘积分期望
重排模型只要求病例的位置和时间,不要求背景人口。它构造期望的方式很巧妙:设 c_zt 是位置 z、时间 t 的病例数,空间边际 C_z = Σ_t c_zt,时间边际 C_t = Σ_z c_zt,总数为 C。如果空间和时间相互独立,期望就是 e_zt = C_z · C_t / C,窗口内的期望就是这些 e_zt 在窗口内的和。这相当于把时间和空间标签在所有病例之间随机打乱后得到的期望——所以叫“重排”。
这个方法的好处是自动调整了纯粹的空間聚集和纯粹的时间聚集:如果某个地区全年病例都多,空间边际已经吸收了这部分,不会被误判成时空交互。坏处也在这里——它没有分母,区分不了“病例多”和“人口多”。说明书里明确提醒:研究期超过一年就要非常小心,如果不同地区人口增长速率不同,重排模型无法区分“因为人多了所以病例多”和“因为风险高了所以病例多”。这个坑在实际分析里非常常见。
3.3 指数、正态与加权正态:连续属性靠置换保住 p 值的有效性
指数模型是为生存时间数据设计的。每个病例有两个时间变量:诊断时间(被扫描的时间维度)和生存时间(病例的一个属性),外加一个 0/1 的删失标记。扫描窗口扫的是诊断时间,检验的是窗口内病例的生存时间是否异常地长或短。似然函数基于指数分布,但 p 值通过置换空间位置和时间标签获得,所以真实生存时间分布不必是指数分布——这是置换检验相对参数检验的核心优势。
正态模型处理的是每个观测带一个连续属性(比如出生体重),检验窗口内该属性的均值是否异常。同一个置换逻辑保证了 p 值对任意连续分布都有效,但有个前提:极端值会强烈影响均值,分析前必须截断或缩尾。加权正态版本用于每个观测方差不同的情况,比如每个区域的平均值基于不同样本量,用样本量或方差倒数做权重。
权重有个性质值得记住:所有权重乘同一个常数,结果不变;所有观测权重相同时,加权版本退化为标准版本。这不是巧合,而是因为似然比只依赖于均值的相对精度,权重的绝对尺度被约掉了。
3.4 有序与多项:把分类变量塞进同一个似然框架
有序模型的类别有顺序(早期、中期、晚期),可以搜“高类别偏多”的聚集,也可以搜“低类别偏多”的聚集,或者两者同时搜。把类别顺序反过来,等价于把高聚集搜索换成低聚集搜索——这个对称性说明有序模型本质上是在用类别的秩信息。
多项模型的类别没有顺序,同时搜索任意类别比例的异常升高或降低。它天然地处理了多重比较:因为类别组合很多,p 值计算里已经考虑了这些组合的数量。实现上两者都是把每个病例当成带位置和类别标签的观测,窗口内类别分布的偏离用多项似然比去量。设计动机很直白:病理亚型、肿瘤分期、教育程度这类变量强行编码成 0/1 会丢信息,用有序或多项模型才能把类别信息完整保留在似然里。
4. 离散与连续的位置观:坐标是已知的,还是被当成随机变量
4.1 离散扫描:位置固定,随机性全在病例计数上
离散扫描统计里,位置来自用户提供的坐标或区域质心,位置本身是非随机变量。零模型只对每个位置上的病例计数做随机假设,位置集合是给定条件。所以离散扫描的随机性在计数上,不在位置上。
这个定义解释了为什么离散扫描能处理汇总数据:每个县、每个邮区是一个位置,坐标可以用质心或人口加权质心。它也解释了为什么离散扫描的结论依赖位置划分的粒度——同一个研究区,按县汇总和按普查区汇总,能检测到的聚集尺度是不一样的。
4.2 连续扫描:点在凸多边形内随机,强度函数才是检验对象
连续扫描的零假设是研究区域内的齐次泊松过程:点可以出现在凸多边形内的任何地方,强度为常数。检验的对象是“是否存在某块子区域,点密度显著高于均匀密度”。窗口是圆或椭圆,面积参与期望计算,μ = C · (窗口面积 / 研究区面积)。
要求研究区是凸多边形集合,是因为窗口需要落在研究区内、面积计算才成立。典型场景是生态学里的鸟巢分布、犯罪学里的案发点分布——位置本身就是研究对象的一部分,不能说“这个位置固定,只是计数随机”。
4.3 背景人口不均与人口插值:把「预期」算对才能谈「偏高」
泊松模型最容易出错的地方在期望值的算法。如果直接用“窗口内病例数占总病例数的比例”来算期望,隐含假设是人口均匀。SaTScan 用每个位置的风险人口来分摊,非均匀背景被吸收进 μ_i 里。
时间维度上还有一个细节:人口数据往往只在普查年份有,研究期却是连续的。SaTScan 做线性插值——相邻普查年之间按时间比例插值,第一个普查年之前用第一年值,最后一个普查年之后用最后一年值。这个处理很实用,但有个后果:如果研究期跨越多个普查年却只给了一个人口数,插值就退化成常数,等价于假设人口不变。这在人口快速增长的地区会系统性高估早期、低估后期的期望病例数,直接影响聚集的判定。
另一个容易被忽略的是人年口径。泊松模型里人口应该按人年算,而不是某个时点的人数。研究期内有人口流动时用中途人口近似会有偏差,尤其是研究期较长的情况。
5. 实践中的坑:p 值跑出来之后,哪些地方最容易翻车
5.1 坑一:跨年数据用重排模型,人口增长被当成疾病增加
现象:跑出来的聚类集中在人口增长快的地区,时间窗口横跨好几年;把背景人口补进去用泊松模型再跑一遍,聚类消失或大幅减弱。
原因:重排模型的期望来自空间边际乘时间边际,没有分母,无法区分“病例多”和“人口多”。说明书里专门提醒研究期超过几年要慎用,原因就在这里——新社区建成带来的病例增加和真实疾病风险上升在数据上长得一模一样。
解决:研究期较长或人口变化明显时改用离散泊松模型并提供人口文件,或者把研究期切成一年以内的片段分别分析,再看跨片段的一致性。
5.2 坑二:正态模型不截断极端值,一个离群点能拖出假聚集
现象:窗口内均值被一两个极端观测拉偏,出现一个相对风险不高但 p 值很小的“聚集”。
原因:正态模型的似然比由窗口内均值驱动,个别离群值会显著改变均值,而置换检验只保证 p 值对分布形状稳健,不保证对极端值稳健。
解决:分析前对连续属性做缩尾或截断,比如把超过 99 分位数的值压到 99 分位;同时检查离群值是不是录入错误。如果每个观测的方差不同,直接上加权正态模型,用样本量或方差倒数做权重。
5.3 坑三:扫描窗口开太大,似然比被少数几个点主导
现象:一个横跨半个研究区、时间上覆盖整个研究期的“聚集”被判显著,但相对风险数值并不高。
原因:最大空间窗口比例和时间窗口设得过大,似然比由大窗口主导,窗口内包含了太多病例,接近把整个研究区吞进去。这时候“聚集”其实是在描述全局,不是局部异常。
解决:根据研究目的设定最大窗口,空间通常不超过总人口的 50%、时间不超过总研究期的一半;同时检查最大窗口附近是否出现边界聚集,如果有,说明窗口约束在起作用,需要放宽或换模型。
5.4 坑四:输入文件对不齐,程序不报错但结果全错
现象:程序正常跑完、没有任何警告,但期望病例数明显不对,或者总病例数对不上手工汇总。
原因:.cas、.pop、.geo 三类文件的位置标识或坐标格式不一致;.pop 的年份、年龄组、性别等分类键与 .cas 不匹配;坐标单位混用(度与米)。这类错误不会触发解析失败,只会让连接悄悄错位。
解决:跑之前先用小样本核对三类文件的键列和分类变量;先看程序输出的总病例数、总人口是否等于手工汇总;坐标统一到同一投影后再导入。
5.5 坑五:把 p 值读成单个窗口的显著性
现象:看到某个聚类的 p=0.03 就下结论“这个区域有聚集”,换个窗口参数后 p=0.06 就改口“没有聚集”。
原因:报告出的是扫描统计量最大值的蒙特卡洛 p 值,它校正的是“所有窗口里最极端的那一个”,不是针对某个固定窗口的单次检验。p 值会随窗口集合、置换次数、数据版本变化,把它当成固定假设检验的 p 值来读会过度解读。
解决:把 p 值和相对风险、窗口内病例数一起读;用不同最大窗口和置换次数做敏感性分析;报告里写清窗口参数和置换次数,让读者能复现。
6. 把显著性做成监测工具:参数选择本身就是一次建模决策
回顾性分析和前瞻性监测的差别,本质上是对“时间”的用法不同。回顾性分析拿到全部时间的数据,可以枚举所有时间窗口;前瞻性监测每来一个新时点只能检查“以当前时点结束”的窗口,而且每一轮都要重算 p 值,整体假阳性率会随检验次数累积。SaTScan 把时间当成扫描的一个维度,天然适合这种逐时点的检查,但参数选择会直接改变结论:窗口上限设得宽,早期就能报警但假阳性多;设得窄,报警晚但更稳。
| 模型 | 必需输入 | 推断方式 |
|---|---|---|
| 离散泊松 | 病例、人口、坐标 | 按人口比例重抽样病例 |
| 伯努利 | 病例、对照、坐标 | 按总体病例比例重抽样状态 |
| 时空重排 | 仅病例、时间 | 置换时间标签 |
| 指数、正态 | 病例、连续属性、坐标 | 置换位置与属性 |
| 连续泊松 | 点坐标、研究区边界 | 在凸多边形内重抽样点 |
另一个取舍是置换次数。999 次只能分辨到 0.001 量级,要报更小的 p 就得加次数,耗时线性增长。如果只关心 p 是否小于 0.05,999 次通常够用;如果要做多重比较校正或报告精确 p 值,就得提高次数。
我自己的血泪经验是:拿到任何一份扫描结果,先别急着看聚类,先手工核对总病例数和总人口,再确认窗口参数和置换次数。有一次我用重排模型跑了三年跨度的数据,聚类全落在新建城区,换成泊松模型加人口文件后聚类消失——那是人口增长,不是疾病上升。那以后我每跑一次时空分析,都先把研究期切到一年以内,切不了就补齐人口文件再用泊松模型,这个习惯再没翻过车。希望帮到你。
本文还有配套的精品资源,点击获取