一位做临床数据的朋友前两天把一张四格表甩给我四个数是 68、32、45、55附一句P 值到底看哪一行。这张表用 SAS 的 PROC FREQ 跑一遍输出里会同时冒出 Pearson Chi-Square、Continuity Adj. Chi-Square、Likelihood Ratio Chi-Square、Mantel-Haenszel Chi-Square 四五个卡方值还有一堆 Fisher 精确检验的结果。新手最常犯的错不是不会写代码——tables a*b / chisq;谁都会敲——而是不知道该读哪一行的 P 值以及什么时候这个 P 值根本不该信。这篇就把 SAS 做卡方检验这件事从头到尾拆一遍列联表怎么造、TABLES 语句有哪些关键选项、2×2 表上的三个大坑、效应量怎么算、分层数据怎么用 CMH 处理、结果怎么用 ODS OUTPUT 批量落到数据集里。内容偏实战代码可以直接抄跑完对照着看输出就能用。1. 从一张四格表说起卡方到底在比什么1.1 用手算把 68/32/45/55 拆开先把这张表摊平。试验组 100 例有效 68、无效 32对照组 100 例有效 45、无效 55。合计 200 例其中有效 113 例、无效 87 例。卡方检验的核心只有一句话如果分组和结局完全无关每个格子里应该有多少人这个应该就是期望频数Expected算法是所在行合计 × 所在列合计 ÷ 总合计。试验组有效格的期望值 100 × 113 ÷ 200 56.5试验组无效格 100 × 87 ÷ 200 43.5对照组两格同样都是 56.5 和 43.5。然后看实际观察值偏离期望值有多远偏离量平方、除以期望值、四个格子加起来格子观察值 O期望值 EO−E(O−E)²/E试验组·有效6856.511.52.3407试验组·无效3243.5−11.53.0402对照组·有效4556.5−11.52.3407对照组·无效5543.511.53.0402合计200200010.762χ² Σ(O−E)²/E 10.762自由度 (2−1)×(2−1) 1。查卡方分布χ² 10.762、df 1 对应的 P 值约 0.001。所以分组与结局无关这个前提被推翻了试验组和对照组的有效率差异有统计学意义。手算这一遍不是为了在项目里炫技而是为了在 SAS 输出的四五个卡方值里能对上号。你后面会看到SAS 打印出来的 Pearson Chi-Square 就是 10.762一个字都不差能对上号就说明你读对了行。对不上号的时候八成是用了加权、插补或者格式重组过的数据那时候要回头查数据而不是查代码。1.2 三条容易被忽略的底线卡方检验的理论基础是大样本下统计量近似服从卡方分布这个近似成立有三个前提缺一个结果就不可靠。期望频数不能太小。常规要求是所有格子的期望频数都 ≥ 5宽松一点是期望频数小于 5 的格子不超过总格子数的 20%。上面那张表期望值都是 43.5 和 56.5放心用。但如果试验组只有 12 个人期望值掉到 3 以下卡方分布的近似就崩了SAS 会在输出里直接给你一句告警。观察值之间要相互独立。一个人只能进一个格子。同一批病人前后测两次塞进四格表、或者把一个家庭的多个成员当成独立个体都会让卡方值虚高、P 值虚低。这是审稿人最喜欢抓的点也是实际项目里最容易犯的错。总样本量要有底线。一般建议总 N 至少 40。N 小于 40 时即使用 Fisher 精确检验检出真实差异的能力也很有限这时候该讨论的是研究设计而不是统计方法。1.3 什么时候不该用卡方有几类情况卡方不是最佳选择而是错误选择。一是配对设计。同一批对象前后测量、或者按年龄性别一一配对的病例对照研究数据是相关的要用 McNemar 检验SAS 里还是 PROC FREQ但要在 TABLES 语句后面的 AGREE 选项上做文章或者用/ AGREE配合配对表结构跟独立样本的卡方完全是两回事。二是结局是有序分类。比如疗效分成痊愈、显效、有效、无效四级直接当无序分类跑卡方会丢掉顺序信息检验效能白白浪费。这种要用趋势检验或者秩和类方法后面第 5 节会讲。三是极小的表。2×2 且期望频数不足直接上 Fisher 精确检验别硬套卡方。四是同一个数据反复做几十次卡方。比如 30 个基因位点各跑一次这时候单个 P 值 0.03 其实说明不了什么需要多重比较校正PROC MULTTEST 或者 Benjamini-Hochberg 法都得安排上。2. 数据进 PROC FREQ 之前列联表到底该怎么造2.1 明细数据和汇总频数表两条路的写法差别实际拿到的数据有两种形态代码写法完全不同混用是新手第一大坑。形态一一行一个观测。200 个病人 200 行每行有trt和outcome两个变量。这时候直接跑就行proc freq datawork.trial_detail; tables trt*outcome / chisq nocol nopercent; run;形态二一行一个格子。数据只有 4 行长这样data work.trial_agg; input trt $ outcome $ cnt; datalines; 试验组 有效 68 试验组 无效 32 对照组 有效 45 对照组 无效 55 ; run;这种汇总数据必须加 WEIGHT 语句告诉 SAS每个格子代表多少个人proc freq datawork.trial_agg; tables trt*outcome / chisq nocol nopercent; weight cnt; run;不加 WEIGHT 会怎样SAS 会把这 4 行当成 4 个人来算N4卡方值小得可怜P 值接近 1结果完全错误但不会报错——这是最阴险的地方。所以每次拿到数据集先跑一句proc freq dataxxx; tables trt*outcome; run;看看 N 和你知道的总样本量对不对。对不上就说明有问题。2.2 权重与缺失值N 对不上的两大元凶第一条元凶就是上面说的 WEIGHT 漏写或者多写。还有个小概率情况数据里既有明细又做了汇总误用了两次权重N 会翻倍。第二条元凶是缺失值。PROC FREQ 默认把任一变量缺失的观测整个排除掉而且不报错只在输出表底部悄悄写一行 Frequency Missing。如果你的结局变量里有 8 个未知状态你的 N 就变成 192 而不是 200各格百分比的分母也跟着变小。项目汇报里分母对不上是硬伤。处理方式有两个方向。要么先明确缺失的归因用where outcome is not missing;显式限定并在报告里说明要么把缺失单独作为一类纳入分析前提是缺失有临床意义、并且要在方法学部分说清楚。我个人的习惯是分析前先跑一遍proc freq; tables _all_ / missing; run;把所有变量的缺失情况一次看全比挨个查要快得多。顺便说 WEIGHT 的一个细节权重可以是小数SAS 照算不误但卡方检验的自由度是基于格子数算的、不是基于权重和算的所以非整数权重只在你明确知道自己在做什么时才用。质量调查里的抽样权重就属于这种用它之前最好确认期望频数是否还站得住。2.3 中文标签、格式与输出整洁度跑出来的表是要给人看的。两个细节能显著提升可读性。第一是变量标签。原始数据里的变量名往往是trt_cd、outcome_fl这种缩写输出表头直接显示变量名非技术同事看不懂。加一句label trt分组 outcome疗效评价;就能解决。列联表的行列表头会直接用标签效果立竿见影。第二是格式。如果trt存的是 1 和 2输出表头会是1和2加一个自定义格式就能变成试验组和对照组proc format; value $trtf 1 试验组 2 对照组; value $outf 1 有效 0 无效; run; data work.trial_detail; set work.trial_detail; format trt $trtf. outcome $outf.; run;格式只在显示层生效不改变底层取值所以后续任何统计过程都不受影响可以放心加。中文还有个常见问题SAS 9.4 在 Windows 上如果是中文版环境数据文件是 UTF-8 而会话编码是 GB18030或者反过来标签会显示成乱码。先用proc options optionencoding; run;看一眼当前会话编码再决定源数据文件怎么存。真遇到乱码也可以先不管把结果导出成 CSV 再转编码但那样中间的核对环节会很难受不如一开始就把编码统一掉。另外提一句环境层面的常识PROC FREQ 属于 Base SAS 的范畴标准部署装完就能用不需要额外模块授权。这也是为什么做描述性统计和列联表分析时它几乎是所有人的第一选择——门槛低、结果稳、语法几十年没大改过。3. TABLES 语句精讲从最小可用到能交付3.1 最简一行代码与它打印出来的每一块内容proc freq datawork.trial_detail; tables trt*outcome / chisq; run;这一跑会打印两大块。第一块是交叉表本身每个格子有四行字Frequency频数、Percent占总数的百分比、Row Pct行百分比、Col Pct列百分比。这里有个很多人栽过的坑汇报有效率应该看 Row Pct因为行是分组、列是结局如果表写成了tables outcome*trt行列表头一换Row Pct 的含义就完全变了68.00 会变成 60.18。这个错误在真实项目里出现过无数次改表方向之前先想清楚谁是分母。第二块是 Statistics for Table of trt by outcome里面列出的统计量包括统计量含义本文例子取值Chi-SquarePearson 卡方就是手算那个10.7622Likelihood Ratio Chi-Square似然比卡方基于对数似然而非平方和10.8497Continuity Adj. Chi-SquareYates 连续性校正卡方仅 2×2 表出现9.8458Mantel-Haenszel Chi-Square基于行秩相关的卡方2×2 时与 Pearson 略有差异10.7084Phi Coefficientφ 系数卡方开方再除以 N0.2320Contingency Coefficient列联系数0.2260Cramers V克拉默 VR×C 表的关联强度0.2320Fishers Exact Test精确检验2×2 表会一并给出P0.0017 左右Pearson 和似然比卡方在大样本下结论通常一致不一致时一般说明某些格子的期望频数偏小需要警惕。3.2 EXPECTED、CELLCHI2、ALPHA 这些选项各自改变什么光有 P 值不能交付。审稿人和老板都会问哪个格子贡献最大期望频数够不够。下面这几个选项就是为回答这些问题准备的。expected让每个格子多打印一行期望频数。这是判断卡方是否适用的最直接依据。输出表里带 Expected 的行一旦有小于 5 的数值就要开始考虑替代方案。deviation打印观察值减期望值的差值正负号直接告诉你哪个格子的实际人数高于/低于预期。cellchi2打印每个格子的 (O−E)²/E也就是卡方值的分项。哪个格子的这个数最大它就是差异的主要来源。上面那个例子里试验组无效格和对照组无效格各贡献 3.04是总卡方 10.76 的主要来源——换句话说差异主要来自无效这一侧而不是有效这一侧。这种解读在写报告时非常有用。alpha0.01只改变置信区间的置信水平不改变 P 值。这是个高频误解很多人以为改了 alpha 就是改了显著性判定标准其实 SAS 只是按这个水平算优势比的置信区间。判定显著性永远是你自己拿 P 值和事先定好的 α 去比。nocol nopercent砍掉列百分比和总百分比只留频数和行百分比。表小一点汇报时不容易看串行。反过来如果分析的是列方向的分母就用norow nopercent。missing把缺失值当成一个有效类别单独列出来前面说过分析前排查数据时常用。noprint抑制交叉表输出只留统计量。批量出结果的时候很省屏幕。3.3 多个变量一起跑TABLES 列表写法与 LIST 布局一个数据集往往有十几个分类变量要两两做卡方逐个写 PROC FREQ 太啰嗦。TABLES 语句支持一次列多张表proc freq datawork.trial_detail; tables (sex agegrp bmi_grp)*outcome / chisq nocol nopercent; run;括号里的变量会分别与 outcome 交叉一次出三张表的统计量。如果还想让它们两两之间也交叉把*改成星号连接的形式或者用指定层级数不过实际项目里很少有人这么干容易把输出撑爆。当表超过两维时比如tables center*trt*outcomeSAS 默认会把每个中心打印成一张独立的 2×2 表同一屏上刷出几十页。这时候加list选项改成清单式布局每行一个格子、带 strata 列看起来清爽很多proc freq datawork.multicenter; tables center*trt*outcome / list nocol nopercent; weight cnt; run;清单布局的另一个好处是容易肉眼比对各中心的分布是否均衡为后面做 CMH 分层分析做铺垫。如果只想看统计量、不想看几十页表就noprint把输出全部导向 ODS 数据集这一招在第 8 节会展开。4. 2×2 表上最容易翻车的三件事4.1 Pearson、Continuity Adj.、Mantel-Haenszel 该选哪个这是被问得最多的问题没有之一。SAS 在 2×2 表的统计量表里同时给出三个卡方它们不是同一个东西Pearson Chi-Square 是不带任何校正的经典卡方也是手算能对上的那一个。Continuity Adj. Chi-Square 是 Yates 连续性校正版算法是把每个格子的偏差先减掉 0.5 再平方χ²_c Σ(|O−E| − 0.5)²/E。上面例子里就是 (11.5−0.5)²/56.5 ×2 (11.5−0.5)²/43.5 ×2 9.846比 10.762 小一圈P 值从 0.001 变成 0.0017。Mantel-Haenszel Chi-Square 走的是另一条路等于 Pearson 值乘以 (N−1)/N在 2×2 表上数值非常接近但不相等。选哪个我的规则很简单总样本量大于 40 且所有期望频数 ≥ 5直接用 Pearson不校正。Yates 校正的思路是让检验更保守代价是真实差异存在时更容易漏掉样本量够的时候没必要自缚手脚。只有当表比较小N 在 20 到 40 之间又不想上精确检验时才把校正版作为折中方案。需要提醒的是有些统计软件在 2×2 表上默认给的是校正后的值如果你用两个工具交叉验证结果看到的卡方值差了零点几先别慌着怀疑代码很可能只是默认设置不同。4.2 期望频数小于 5 的告警怎么处理SAS 触发告警的条件是期望频数小于 5 的格子所占比例超过 20%输出里会出现类似这样一行WARNING: 25% of the cells have expected counts less than 5. Chi-Square may not be a valid test.看到告警先做三件事按顺序来。第一步确认是哪个格子出了问题。加expected cellchi2跑一遍把带 Expected 的那一行扫一遍。很多时候只是某一格特别薄其他格子都很健康。第二步判断能不能合并类别。如果某个格子是仅有 2 例有效的极稀有类别而且和相邻类别的临床意义接近合并后问题就没了。合并要有学科依据不能纯粹为了凑期望频数硬合。第三步不能合并就换方法。总样本量还过得去、只是分布偏斜的可以用 Fisher 精确检验分层很细、每层例数都很少的考虑用精确 logistic 回归。千万不要为了绕开告警去人为加几例数据那是数据造假。顺带一个经验连续型变量分组后做卡方如果分组是按中位数二分两组例数会很接近期望频数通常没问题但如果按临床界值分很容易一边倒比如按 5.6 的界值切出来一组只有 9 例卡方就悬了。分组之前先跑一句proc freq; tables 分组变量; run;看看分布能省掉后面很多返工。4.3 Fisher 精确检验什么时候必须上Fisher 精确检验不依赖大样本近似直接基于超几何分布计算当前这张表或更极端情况出现的精确概率所以小样本下它是标准答案。什么时候上总 N 小于 40或者 N 大于 40 但有格子期望频数小于 5或者有格子观察频数为 0。另外在 2×2 表上即使你只写了/ chisqSAS 一般也会同时打印 Fisher 的结果不同版本、不同设置下如果没自动出就在后面补一句exact fisher;显式指定所以判断成本几乎为零看 P 值哪个更保守、更符合方法学要求就行。代价要说清楚三点。第一精确检验是条件检验结论只对当前这一张表的行合计和列合计成立外推性比卡方弱。第二表一大就慢3×3 以上、格子又多的时候计算量指数级上升SAS 会报内存或时间超限需要设exact fisher / maxtime120;之类的参数控制。第三真实项目里如果 N 只有 30 例还要比两组有效率我一般会劝合作者把结论写成提示性结果而不是硬下统计学判断——检验方法再精确也补不回信息量的缺失。5. P 值之外效应量与有序趋势检验5.1 RELRISK、RISKDIFF、ODDSRATIO 与方向陷阱P 值告诉你有没有差异效应量告诉你差多少。这在 2×2 表上是刚需因为 N 一大0.5% 的差异也能做出显著。加几个选项就能拿到全套proc freq datawork.trial_detail; tables trt*outcome / chisq relrisk riskdiff oddsratio nocol nopercent; run;relrisk给出相对危险度本例 0.68/0.45 1.511和它的置信区间另有列 1/列 2 两个方向的两套值。riskdiff给出风险差0.68−0.45 0.23及其区间。oddsratio给出优势比(68×55)/(32×45) 2.597含义是试验组有效的优势是对照组的 2.6 倍。方向陷阱是这一块的高频事故SAS 默认按行变量的第一水平比第二水平、列变量的第一水平比第二水平来算。如果你的行变量是试验组/对照组、列变量是有效/无效那 OR 表达的是有效的优势比一旦数据里 0/1 编码反了或者用 format 做了反向映射OR 就变成倒数变成 0.385从有效 2.6 倍变成无效 2.6 倍含义完全反了。核对方法很简单拿手算值对一遍对不上就是方向反了。5.2 MEASURES 选项给出的一整套效应量不确定该报哪个效应量的时候measures是个偷懒但好用的选择它一次性打印一大张表包含 Pearson 相关系数、Spearman 相关系数、Lambda、不确定系数、Gamma、Somers D、Kendalls tau-b、Stuarts tau-c 等等。这张表信息量很大但直接往论文里搬是错的。挑选原则看变量性质两个无序分类变量报 Cramers V 或 φ两个有序分类变量报 Kendalls tau-b 或 Gamma一个有序一个无序Somers D 比较合适。有序变量用无序的关联指标会低估关联强度白白丢掉信息。还有一个细节这些系数各有取值范围和解释惯例。φ 和 Cramers V 在 0 到 1 之间0.2 左右属于弱关联Gamma 和 Kendalls tau-b 在 −1 到 1 之间带方向。汇报时最好同时给出系数和它的解读标准不然读者没法判断 0.23 算强还是算弱。5.3 有序分组的趋势检验Cochran-Armitage如果列联表的一维是有序的比如剂量组分成低、中、高或者年龄分成 40、40-60、60那么是否存在随剂量递增的线性趋势比三组之间有没有差异更有价值。SAS 里的写法是加trend选项proc freq datawork.dose_trial; tables dosegrp*response / trend nocol nopercent; weight cnt; run;前提是其中一维是二分类比如应答/未应答另一维是有序多分类。SAS 会输出 Cochran-Armitage 趋势检验的 Z 统计量和单侧、双侧 P 值。这里有个关键点容易被忽略趋势检验的得分赋值方式会直接影响结果。剂量是 10mg、20mg、50mg 这种不等距的等距赋分1、2、3就等于假设间隔相同会稀释真实的剂量效应。不等距时可以在 TABLES 语句里用 SCORES 选项指定或者把实际的剂量值作为一个数值变量传给 SAS 用。我个人在做剂量-反应分析时习惯用实际剂量值而不是编组号这样得到的趋势检验更有解释力。另外别忘了趋势检验和普通卡方不是替代关系而是互补关系。两种都报出来更稳妥卡方不显著但趋势检验显著说明存在方向性效应但组间两两比较没有差异卡方显著但趋势检验不显著说明差异不是单调递增的可能是某个中间剂量组异常。6. 分层与匹配CMH 检验解决的是混杂问题6.1 为什么把各中心的数据直接合并会出问题多中心研究是最典型的场景。三个中心各自做出来的有效率可能差距很大原因可能在于中心之间病人基线不一样——A 中心收的都是轻症C 中心重症偏多。这种情况下把三个中心的数据简单相加再做卡方会引入中心这个混杂因素得到的关联估计是有偏的。更麻烦的是辛普森悖论完全可能出现每个中心内部都是试验组更好合并之后反而对照组更好的情况。我确实在真实数据上见过一次原因是各中心的样本量差异极大中心本身的效应被样本量加权后盖过了处理效应。所以只要数据里有分层变量中心、性别、年龄段、疾病亚型就应该考虑先分层分析再看能不能合并。6.2 CMH 三种统计量分别在检验什么Cochran-Mantel-Haenszel 检验的做法是在每一层内部计算这层的关联再按权重把所有层的结果合并成一个总的检验统计量。SAS 的写法就是在原有的表结构前面加一层分层变量proc freq datawork.multicenter; tables center*trt*outcome / cmh nocol nopercent; weight cnt; run;输出里的 CMH 统计量有三行对应三种不同的备择假设统计量检验的是适用场景Correlation行、列变量都是有序时的线性相关2×2 分层表最常用等价于合并优势比检验Row Mean Scores Differ行变量有序、列变量无序剂量分组、疾病分期与分类结局General Association行列都无序最通用的关联检验任何分层表都能用2×2×K 的分层表上一般直接看 Correlation 那一行它的 P 值就是控制分层因素后处理与结局是否有关联。在同一段输出里SAS 通常还会一并给出合并的 Mantel-Haenszel 优势比及置信区间这比 OR 本身更有汇报价值因为它是控制混杂后的效应量不同的版本和选项下这些结果的摆放位置会略有差别找不到就往下翻几屏。6.3 Breslow-Day 齐性检验不通过时的处理合并的前提是各层的效应方向一致、大小接近这个前提用 Breslow-Day 齐性检验来验证。P 值大于 0.05 说明各层同质可以放心合并报一个合并 OR 就够。P 值小于 0.05 说明层间效应不一致这时候报合并 OR 就是把不同东西平均了没有意义。处理思路有三条。第一做交互检验把分层变量和暴露变量的交互项放进回归模型明确指出处理效应在不同中心之间不同这本身就是有价值的发现。第二如果目标只是控制混杂而不是追求效应量可以改用 Mantel-Haenszel 方法给出调整后的检验结果同时说明层间异质性。第三只有两三层且例数充足时逐层报告 OR 和置信区间把它们画成森林图读者一眼就能看出哪一层的效应偏离了。我一般会同时报齐性检验和逐层结果即便同质也保留分层表这是审稿时最容易被追问的地方提前备好省得来回补。顺便分层变量的水平数别太多每层例数太薄的时候MH 估计和齐性检验都不稳通常建议每层至少 10 到 20 例。7. 拟合优度与精确检验单向表的写法和用法7.1 单向表 TESTP 做拟合优度卡方不只能比两组。单向表只有一个变量上的卡方做的是拟合优度检验观察到的各类别频数和某个理论分布或期望比例是否一致。data work.survey; input answer $ cnt; datalines; 非常满意 86 满意 124 一般 62 不满意 28 ; run; proc freq datawork.survey; tables answer / chisq testp(40 30 20 10); weight cnt; run;TESTP 后面的比例按类别顺序一一对应这里表示理论上期望的比例是 40%、30%、20%、10%。注意顺序TESTP 的取值是用 TABLES 语句里类别出现的顺序来对应的如果类别不是按你预期的字母或数值顺序排列比例就会被套错。稳妥做法是先用tables answer / list;看一眼各类别的实际顺序再写 TESTP。还有个小细节SAS 要求这些比例加起来等于 1用百分数就是 100差一点点它可能有提示所以写 1/3 这种循环小数的时候用 (16.667 16.667 16.667 16.667 16.667 16.665) 这种凑数写法最省事。如果只是想检验各类别是否均匀分布连 TESTP 都不用写默认就是等比例假设。7.2 EXACT 语句的常见组合与开销EXACT 语句写在 PROC FREQ 里、TABLES 之外常见组合有几种proc freq datawork.small_trial; tables trt*outcome / chisq nocol nopercent; exact fisher chisq; weight cnt; run;exact fisher只要 Fisher 的精确 P 值exact chisq会给出卡方统计量的精确分布下的 P 值。还有exact or可以给出优势比的精确置信区间小样本研究报告里非常吃香因为 Wald 型区间在小样本下偏差明显。性能上要有个概念精确计算的开销随表维度增长得非常快。2×2 表几乎瞬间完成3×3 表还算能接受4×4 以上或者分层表就很容易跑到超时。真需要跑大表精确检验加maxtime控制单次计算时间或者先用exact chisq / mc;用蒙特卡洛模拟估算 P 值——虽然结果带一点随机波动设好种子就能复现实际项目里完全够用。7.3 单向表和双向表在语法上的差别单向表和双向表在语法层面的差别不大但很关键对比项单向表双向表TABLES 写法tables var;tables var1*var2;CHISQ 的含义拟合优度检验独立性检验TESTP 选项可用指定期望比例不适用Fisher 精确检验不适用2×2 常用输出中的 ODS 表名OneWayChiSqChiSq / FishersExact单向表的 CHISQ 输出里统计量叫Chi-Square Test for Specified Proportions指定了 TESTP 时或Chi-Square Test for Equal Proportions没指定时这个标题一定要看清楚它直接说明了检验的零假设是什么。做过一次项目就发现有人拿着Equal Proportions的结果去汇报符合预设比例差得挺远。8. 把结果搬进数据集ODS OUTPUT 与宏循环8.1 用 ods trace 找到表名再用 ODS OUTPUT 抓数前面所有操作都在看屏幕输出但一个项目几十上百张表手抄进 Excel 是不现实的。正确做法是把统计量直接落到数据集里。首先要找到 ODS 表名。SAS 的每个过程在 ODS 系统里都注册了若干张表用ods trace就能看到ods trace on; proc freq datawork.trial_detail; tables trt*outcome / chisq relrisk oddsratio; run; ods trace off;日志里会打印出每一张输出表的名称卡方统计量的表名通常是 ChiSqFisher 的是 FishersExact效应量的在 Measures 里CMH 的在 CMH 里单向表拟合优度的在 OneWayChiSq 里交叉表的原始格子数据在 CrossTabFreqs 里。表名在不同版本之间基本稳定但万一抓不到跑一次 trace 是最快的确认方式。拿到表名就能定向抓取ods output ChiSq work.chisq_out Measures work.measures_out; proc freq datawork.trial_detail; tables trt*outcome / chisq relrisk oddsratio nocol nopercent; run; ods output close;ChiSq 数据集的结构很简单Statistic、DF、Value、Prob 四列每行一个统计量。取值时可以直接where Statistic Chi-Square;。注意ods output close;一定要写否则影响会一直挂着后面所有 PROC FREQ 的输出都会被重定向排查起来很费时间。8.2 宏循环批量跑几十个变量手头有 30 个基因位点要做卡方一个一个写不现实。用宏循环包一层把每个变量的结果分别存到数据集里再合并%macro batch_chisq(dsn, ylist, x); %local i v; %do i 1 %to %sysfunc(countw(ylist, %str( ))); %let v %scan(ylist, i, %str( )); ods output ChiSq work.cs_v; proc freq datadsn; tables x*v / chisq nocol nopercent; run; ods output close; %end; %mend batch_chisq; %batch_chisq(dsnwork.trial_detail, ylistsite1 site2 site3 site4 site5, xtrt);跑完之后用集合引用把散落的数据集拼起来并且用indsname把变量名带上data work.all_chisq; set work.cs_: indsnamedsn; length var $32; var scan(dsn, 2, .); if Statistic Chi-Square; run; proc print datawork.all_chisq noobs; var var Value Prob; run;这一步做完几十个位点的卡方值和 P 值就在一张表里了可以直接where Prob 0.05筛显著项或者用 PROC MULTTEST 做 FDR 校正。手工做这一步的话一个下午就过去了。8.3 结果核对与几个高频问题批量跑完之后必须核对不能直接交。三个核对动作第一抽样复查。随机挑三五个变量手工跑一次 PROC FREQ比对卡方值和 P 值是否一致。数量级一样但小数位差一点通常是抓错了统计量行差得离谱就是变量标签或者权重出了问题。第二检查数据集条数。数据集的观测数应该等于变量数 × 每个变量对应的统计量行数。少了说明某个 PROC FREQ 中途报错日志里搜一下 ERROR 和 WARNING 就能找到。第三检查 P 值的分布。如果几十个变量跑出来 P 值全都集中在 0.9 附近那肯定有问题——要么变量和分组完全独立也就是数据被人为构造过要么两组标签搞反了。真实的生物学或临床数据不会这么整齐。最后说一个高频报错ERROR: Variable xxx not found。宏循环里最常见的原因是变量名拼写和数据集里不一致或者某个变量全是缺失被 SAS 处理掉了。先用proc contents dataxxx; run;把变量清单拉出来对一遍比逐行查代码快得多。另一个常见的是WARNING: Apparent symbolic reference not resolved一般是宏变量名写错或者作用域没传进来在宏里加一句%put v;打出来看看就清楚了。卡方检验本身不复杂难的是每次都问自己三个问题这张表是怎么造出来的、期望频数够不够、输出的那几行到底该读哪一行。这三问答清楚了SAS 里剩下的都是几行代码的事。