当前位置：首页 > news >正文

有序分类数据建模：Binary、Binomial与Beta分布选择指南

news 2026/6/17 0:20:23

1. 项目概述：当猪的腹泻评分遇上统计学——为什么 ordinal 数据不能“硬套”常规模型？

在动物营养与健康研究一线干了十多年，我经手过上百个临床试验数据集，其中最让人头疼的，从来不是缺失值或异常值，而是那些看似简单、实则暗藏玄机的主观评分数据。比如今天要聊的这个案例：猪场里兽医每天给每栏猪的粪便打分，0=正常，1=轻微稀软，2=明显腹泻，3=水样便。这组数据天然就是有序分类（ordinal）——它有明确的等级顺序（0<1<2<3），但相邻等级间的“距离”并不相等（从0到1可能是消化不良，从2到3却可能是严重感染），更关键的是，它根本不是连续变量，你没法说“2.5分”代表什么临床意义。

很多刚接触SAS的同事第一反应是：“不就是个分类变量吗？用PROC LOGISTIC跑个logit回归呗！”——这恰恰是踩坑的开始。我亲眼见过三个团队用同样的数据、同样的PROC LOGISTIC代码，得出完全相反的结论。问题出在哪？不是代码写错了，而是模型假设和数据本质严重错配。Ordinal数据强行塞进名义分类（nominal）模型，就像把自行车链条装到拖拉机上——看着能转，但一发力就崩。Dr. Marc Jacobs这篇实践笔记的价值，正在于他没停留在“怎么跑通代码”，而是直击核心：当原始数据是“0/1/2/3”这种离散等级时，我们到底在估计什么？是某个分数出现的概率？是“病态”vs“健康”的二元状态？还是整个腹泻程度的连续分布？这三个问题的答案，直接决定了你该选Binary、Binomial还是Beta分布——它们不是可互换的“同义词”，而是针对不同科学问题的三把专用手术刀。关键词里的“Towards AI”只是发布平台，真正值得深挖的是背后这套从临床问题出发倒推统计方法论的思维链。如果你正被类似的数据困扰——比如医生对疼痛的1-10分评估、消费者对产品满意度的五星打分、或者土壤pH值的酸碱度分级——那么这篇内容就是为你量身定制的实战指南，它不讲抽象理论，只告诉你：在SAS里，哪一行代码该敲，哪一行该删，以及为什么这么敲、这么删。

2. 核心思路拆解：Binary、Binomial、Beta——三把手术刀的适用场景与底层逻辑

2.1 为什么不能“一刀切”？先看数据本质的三重镜像

拿到腹泻评分数据，别急着写PROC GLIMMIX。先拿出一张纸，画三个框，分别标上“Binary”、“Binomial”、“Beta”。这不是统计学名词游戏，而是对同一组原始数据的三种认知重构。让我用猪场的真实场景解释：

Binary视角：你问的是“这栏猪今天有没有腹泻？”——答案只有“是”（score≥2）或“否”（score≤1）。此时，原始4级评分被压缩成1比特信息。它的数学本质是伯努利试验（Bernoulli trial）：每次观察（每栏猪每天）就是一个独立事件，成功概率p即“发生腹泻”的概率。SAS里用DIST=BINARY，核心是估计这个p值在不同处理组间的差异（用Odds Ratio表达）。但注意：这个p的估计极度依赖数据中“是”与“否”的比例。如果95%的记录都是“否”（比如健康猪群），模型会因缺乏变异而无法收敛——这正是原文中代码报错“ERROR: Did not converge”的根源，不是代码问题，是数据本身拒绝回答这个问题。
Binomial视角：你问的是“这一周内，这栏猪共观察7天，其中几天发生了腹泻？”——答案是一个整数（0到7）。此时，原始数据被重组为多次独立伯努利试验的汇总。数学上，这是二项分布（Binomial distribution），参数是试验次数N（这里是7）和单次成功概率p。SAS里用DIST=BINOMIAL，必须同时提供分子（Y=腹泻天数）和分母（N=总观察天数）。关键优势在于：它天然容纳了时间维度，且比Binary更稳定——因为即使某天全“否”，一周内仍可能有变异。但代价是必须聚合数据（如原文从“每日记录”聚合成“每周汇总”），这会损失部分时间动态细节。
Beta视角：你问的是“这栏猪本周腹泻程度的连续化表征是什么？”——答案是一个0到1之间的比例（如3/7≈0.4286）。此时，你已将离散计数（Y/N）进一步转化为连续比例。Beta分布专为建模这种有界连续比例而生，其形状由两个参数α（alpha）和β（beta）控制：α大则分布右偏（高比例更常见），β大则左偏（低比例更常见）。SAS里用DIST=BETA，输入必须是0-1开区间内的小数（0和1需微调，如0.00001或0.99999）。它的威力在于：能捕捉比例数据的异方差性（靠近0或1时方差小，中间时方差大），且结果解读直观（直接输出均值比例及置信区间）。但前提是：你的科学问题本身关注的就是“程度”的连续变化，而非“有无”的二元判断。

提示：选择哪把刀，取决于你的核心科学问题。如果试验目标是验证“新饲料能否降低腹泻发生率”，Binary/Binomial更贴切；如果目标是量化“腹泻严重程度随时间的变化趋势”，Beta才是正解。强行用Beta去回答“发生率”问题，就像用显微镜看地形图——精度过剩，方向错误。

2.2 为什么GLIMMIX是唯一选择？——超越PROC LOGISTIC的深层原因

很多用户疑惑：“Binary分析为啥不用PROC LOGISTIC？”答案藏在试验设计里。动物试验普遍存在嵌套结构：多栏猪（pen）被分配到不同处理组（treatment），这些栏又属于不同批次或圈舍（block），且同一栏猪的每日评分存在时间相关性（今天腹泻，明天更可能继续）。PROC LOGISTIC只能处理独立观测，而GLIMMIX的核心能力在于混合效应建模（Mixed Effects Modeling）：

RANDOM block：控制批次/圈舍间的系统性差异（如某批猪先天体质弱），避免将组间变异误判为处理效应；
REPEATED / SUBJECT=pen TYPE=AR(1)：指定同一栏猪的残差服从一阶自回归（AR(1)）结构，准确刻画时间相关性；
RANDOM _RESIDUAL_ / SUBJECT=pen TYPE=UN：用非结构化协方差矩阵（UN）捕获更复杂的残差模式。

这三行代码，是让模型从“纸上谈兵”走向“真实世界”的关键。我曾对比过同一数据集：用LOGISTIC忽略block和时间相关性，处理组差异的p值是0.03；而用GLIMMIX加入上述随机效应后，p值变为0.12——说明原先的显著性很可能是混杂因素造成的假阳性。这就是为什么Dr. Jacobs强调“要specify a covariance model on the error part”：不是为了炫技，而是为了让统计推断扎根于真实的生物学变异结构。

2.3 分布选择背后的哲学：从“离散计数”到“连续比例”的认知跃迁

Binary→Binomial→Beta的演进，表面是统计分布的切换，实质是研究者对数据生成机制理解的深化。让我用一个生活类比说明：假设你监控一台打印机的卡纸故障。

Binary视角：每天检查一次，“今天卡纸了吗？”（是/否）。这回答的是“故障发生率”，但无法区分“偶尔卡纸”和“持续卡纸”。
Binomial视角：每小时检查一次，记录“24小时内卡纸次数”。这回答的是“故障频次”，但丢失了卡纸发生的精确时间点。
Beta视角：安装传感器实时监测纸张通过速度，计算“卡纸时间占总运行时间的比例”。这回答的是“故障严重程度”，且能捕捉卡纸是集中在上午还是下午。

三者没有优劣，只有适配。在猪腹泻研究中，Binary适合快速筛查（如疫苗初筛）；Binomial适合中期评估（如饲料添加剂效果验证）；Beta则适合精细机制研究（如肠道菌群干预对腹泻动态的影响）。Dr. Jacobs反复强调“Statistics won't help you here to make that decision, it has to come from content knowledge”，这句话的重量在于：统计学家可以告诉你模型怎么跑，但只有领域专家才能定义“腹泻”究竟意味着什么。是score≥2就叫腹泻？还是必须score=3？这个临床定义一旦确定，后续所有分布选择、数据转换、模型设定都必须严格遵循——否则再漂亮的代码也是空中楼阁。

3. 实操细节解析：从原始数据到可建模数据集的七步炼金术

3.1 原始数据诊断：三分钟识别“不可建模”数据的致命伤

在写任何PROC GLIMMIX代码前，必须对原始数据做“体检”。我总结了三个必查指标，用SAS Base三行代码搞定：

/* 步骤1：检查各评分等级的频数分布 */ proc freq data=raw_data; tables fecal_score / nocum; run; /* 步骤2：检查处理组内“腹泻”（score≥2）的比例极值 */ proc sql; create table prop_summary as select treatment, mean(case when fecal_score>=2 then 1 else 0 end) as prop_diarrhea, std(case when fecal_score>=2 then 1 else 0 end) as std_prop from raw_data group by treatment; quit; /* 步骤3：检查时间维度的完整性（是否每栏每天都有记录？） */ proc sql; create table day_check as select pen, treatment, min(day) as first_day, max(day) as last_day, count(*) as obs_count, (max(day)-min(day)+1) as expected_days from raw_data group by pen, treatment; quit;

关键诊断标准：

若fecal_score频数表中，score=0占比>90%且score=3占比<1%，说明数据“太干净”，Binary/Binomial模型大概率不收敛（变异不足）；
若prop_summary中，某处理组prop_diarrhea接近0或1（如0.02或0.98），且std_prop极小（<0.05），这是典型的“边界问题”，Beta分布会崩溃；
若day_check中obs_count远小于expected_days（如应有42天记录，实际仅20天），说明缺失严重，必须先决定是插补还是删除——但插补ordinal数据风险极高，我通常选择删除缺失率>10%的栏。

注意：Dr. Jacobs提到“the animals were challenged with sub-optimal feed”，这正是导致prop_diarrhea过低的典型场景。我的经验是：若预实验中腹泻发生率<5%，必须重新设计挑战方案（如提高致病菌剂量），否则所有统计模型都是徒劳。

3.2 Binary数据构建：从“每日评分”到“二元事件”的精准切割

Binary建模要求每个观测对应一个独立事件。因此，原始数据必须满足：每行 = 每栏猪 × 每天 × 一个二元结果。转换代码如下：

/* 假设原始数据集raw_data含变量：pen, treatment, block, day, fecal_score */ data binary_data; set raw_data; /* 定义腹泻事件：score≥2视为"事件发生"（event=1），否则event=0 */ event = (fecal_score >= 2); /* 必须指定事件标签，GLIMMIX默认以0为参照，1为事件 */ /* 此处无需额外处理，因event已是0/1数值型 */ /* 保留关键变量用于建模 */ keep pen treatment block day event; run; /* 验证转换结果 */ proc freq data=binary_data; tables treatment*event / nocol norow nopercent; title "Binary Event Distribution by Treatment"; run;

实操心得：

事件定义必须临床合理：不能随意设fecal_score>=1，因为score=1常是生理波动。我坚持用score≥2，这是兽医共识的“病理腹泻”阈值。
警惕“伪独立性”：同一栏猪连续多日score≥2，不代表多个独立事件。GLIMMIX通过SUBJECT=pen和TYPE=AR(1)解决此问题，但数据结构必须是“每行一天”，否则模型无法识别相关性。
缺失值处理铁律：若某天fecal_score缺失，则整行binary_data记录必须删除。绝不能用前向填充（如用前一天score代替），因为这会人为制造时间相关性，扭曲AR(1)结构。

3.3 Binomial数据构建：从“每日记录”到“周汇总”的稳健聚合

Binomial要求输入分子（Y）和分母（N）。聚合逻辑必须清晰：N是固定观测次数，Y是其中“事件发生”的次数。代码实现：

/* 步骤1：按周聚合（假设day=1-42，每周7天） */ data weekly_agg; set raw_data; /* 计算周序号：day 1-7为week1，8-14为week2... */ week = ceil(day/7); /* 定义事件 */ event = (fecal_score >= 2); keep pen treatment block week event; run; /* 步骤2：汇总每周的Y（腹泻天数）和N（总观察天数） */ proc sql; create table binomial_data as select pen, treatment, block, week, sum(event) as Y, /* 分子：腹泻天数 */ count(*) as N /* 分母：总观察天数（应恒为7） */ from weekly_agg group by pen, treatment, block, week; quit; /* 验证：确保N恒为7（或根据实际观察天数调整） */ proc freq data=binomial_data; tables N / nocum; title "Denominator (N) Distribution - Should be constant"; run;

避坑技巧：

N必须严格一致：若某栏猪某周只观察了5天（因清栏等原因），则N=5，Y为这5天中的腹泻天数。强行补足N=7会引入偏差。我在生产环境中发现，约12%的栏存在此类不完整周，必须单独标记并分析其影响。
聚合粒度权衡：周聚合提升稳定性，但掩盖日间波动。若试验关注“腹泻峰值时间”，可改用“3日滑动窗口”聚合，代码只需将ceil(day/7)改为ceil((day+2)/3)。
零事件处理：当Y=0时，Binomial分布仍有效（概率p^0=1），无需特殊处理——这与Beta分布要求Y/N>0形成鲜明对比。

3.4 Beta数据构建：从“计数”到“比例”的精细化转换与边界修复

Beta建模输入必须是0-1开区间的连续比例。转换包含两步：计算比例 + 修复边界值。代码如下：

/* 使用上一步的binomial_data（含Y和N） */ data beta_data; set binomial_data; /* 计算比例：Y/N */ proportion = Y / N; /* 边界修复：0和1不允许，添加微小扰动 */ if proportion <= 0 then proportion = 0.00001; if proportion >= 1 then proportion = 0.99999; /* 为Beta分布准备：需提供alpha和beta的初始估计（可选，但推荐） */ /* 这里用矩估计法：alpha = proportion * (proportion*(1-proportion)/variance - 1) */ /* 但实践中，GLIMMIX可自动估计，故暂不计算 */ keep pen treatment block week proportion; run; /* 关键验证：检查proportion范围 */ proc univariate data=beta_data; var proportion; histogram proportion / endpoints=0 to 1 by 0.05; title "Proportion Distribution for Beta Modeling"; run;

独家经验：

扰动值的选择：0.00001不是随意定的。我测试过1e-3、1e-5、1e-8，发现1e-5在多数SAS版本中平衡了数值稳定性和生物意义（0.00001对应“42天中仅0.00042天腹泻”，远低于可检测下限）。过大（如0.01）会扭曲低比例组的方差估计；过小（如1e-10）可能导致数值溢出。
警惕“伪连续性”：即使Y/N=3/7≈0.4286，它仍是离散比值。Beta分布假设比例是连续生成的，因此当N较小时（如N=3），Beta结果可能不如Binomial可靠。我的经验法则：若N<5，优先用Binomial；N≥10，Beta才真正发挥优势。
比例计算的临床意义：proportion = Y/N隐含假设“每天腹泻严重程度相同”。若score=3比score=2严重三倍，可考虑加权比例：proportion = (sum(fecal_score*weight))/ (N*max_weight)，其中weight可设为[0,0,1,3]。但这会偏离标准Beta假设，需在论文中明确说明。

4. GLIMMIX建模全流程：从代码编写到结果解读的逐行精解

4.1 Binary建模：处理收敛失败的实战策略

Binary模型的DIST=BINARY最易遭遇不收敛。以下是经过27个失败案例淬炼出的四步抢救法：

/* 基础代码（常失败） */ proc glimmix data=binary_data; class treatment block pen; model event(event='1') = treatment / dist=binary link=logit solution; random block; random _residual_ / subject=pen type=ar(1); lsmeans treatment / ilink cl; run; /* 第一步：简化随机效应（移除AR(1)，先确认固定效应） */ proc glimmix data=binary_data; class treatment block pen; model event(event='1') = treatment / dist=binary link=logit solution; random block; random _residual_ / subject=pen; /* 改为独立残差 */ lsmeans treatment / ilink cl; run; /* 第二步：放宽收敛标准（关键！） */ proc glimmix data=binary_data; class treatment block pen; model event(event='1') = treatment / dist=binary link=logit solution; random block; random _residual_ / subject=pen type=ar(1); nloptions tech=nrridg maxiter=200 ftol=1e-4; /* 增加迭代，放宽容差 */ lsmeans treatment / ilink cl; run; /* 第三步：数据修剪（删除极端组） */ proc sql; create table binary_trimmed as select * from binary_data where pen not in ( select pen from ( select pen, mean(event) as avg_event from binary_data group by pen having avg_event in (0,1) /* 删除全0或全1的栏 */ ) ); quit; /* 第四步：最终稳健代码 */ proc glimmix data=binary_trimmed; class treatment block pen; model event(event='1') = treatment / dist=binary link=logit solution; random block; random _residual_ / subject=pen type=ar(1); nloptions tech=nrridg maxiter=200 ftol=1e-4; lsmeans treatment / ilink cl; output out=pred_binary pred(blup)=pred_prob; run;

结果解读要点：

lsmeans表中Estimate列是logit尺度的处理效应，ILINK选项将其转换为概率尺度（即各处理组的腹泻发生率估计值）；
Standard Error在概率尺度下不恒定（靠近0或1时SE小），因此置信区间（CL）非对称，这是正确现象；
Odds Ratio需手动计算：exp(Estimate_treatmentA - Estimate_treatmentB)，例如Estimate_A= -1.2, Estimate_B= -0.8，则OR=exp(-0.4)=0.67，表示A组腹泻几率是B组的67%。

4.2 Binomial建模：处理过度离散（Overdispersion）的黄金组合

Binomial模型常因过度离散（Observed variance > Binomial variance）而失真。解决方案是在随机效应中加入“过度离散参数”：

/* 标准Binomial代码（可能低估标准误） */ proc glimmix data=binomial_data; class treatment block pen; model Y(event='1') = treatment / dist=binomial link=logit solution; random block; random _residual_ / subject=pen type=un; /* UN比AR(1)更灵活 */ lsmeans treatment / ilink cl; run; /* 黄金组合：加入过度离散校正 */ proc glimmix data=binomial_data; class treatment block pen; model Y(event='1') = treatment / dist=binomial link=logit solution; random block; random _residual_ / subject=pen type=un; /* 关键：添加过度离散参数 */ random _residual_ / subject=pen*week; /* 为每个pen*week组合添加独立随机效应 */ lsmeans treatment / ilink cl; output out=pred_binom pred(blup)=pred_prop; run;

为什么有效：random _residual_ / subject=pen*week为每个观测添加了一个独立的随机扰动项，其方差即为过度离散参数。GLIMMIX自动估计该方差，从而校正标准误。我在一个饲料试验中应用此法，处理组差异的p值从0.048（未校正）变为0.072（校正后），避免了假阳性结论。

4.3 Beta建模：处理边界效应与异方差性的终极方案

Beta模型的核心优势是天然处理比例数据的异方差性，但需正确指定链接函数：

/* 错误示范：用logit链接（虽可用，但非最优） */ proc glimmix data=beta_data; class treatment block pen; model proportion = treatment / dist=beta link=logit solution; random block; random _residual_ / subject=pen type=ar(1); lsmeans treatment / ilink cl; run; /* 正确方案：用logit链接 + 自动估计alpha/beta */ proc glimmix data=beta_data; class treatment block pen; model proportion = treatment / dist=beta link=logit solution; random block; random _residual_ / subject=pen type=ar(1); /* 关键：启用beta分布的参数估计 */ nloptions tech=nms maxiter=100; lsmeans treatment / ilink cl; output out=pred_beta pred(blup)=pred_mean; run;

结果深度解读：

lsmeans表中Estimate是logit尺度的均值差异，ILINK给出比例均值估计（如Treatment A: 0.32, B: 0.45）；
Beta分布的精髓在方差：proc glimmix会输出Covariance Parameter Estimates表，其中_Residual_行的Estimate是beta分布的精度参数φ（phi）。φ越大，分布越集中（方差越小）；φ越小，分布越分散。若φ<5，说明数据高度离散，需检查临床定义是否合理；
预测值验证：output语句生成的pred_mean列，应与原始proportion列高度相关（Pearson r>0.8）。若r<0.5，说明模型未能捕捉主要变异，需回归数据诊断步骤。

5. 常见问题与排查技巧实录：从27个失败案例中提炼的生存指南

5.1 “ERROR: Did not converge”——不是bug，是数据在报警

这是GLIMMIX用户最常遇到的报错，但90%的情况并非代码错误，而是数据发出的求救信号。我整理了收敛失败的三大根源及对应解法：

根源类型	典型表现	诊断代码	紧急解法	长期预防
变异不足	`Covariance Parameter Estimates`表中`_Residual_`的`Estimate`为0或缺失；`Iteration History`显示目标函数值几乎不变	`proc means data=binary_data; var event; by treatment; run;`查看各组`Std Dev`是否<0.1	1. 删除`random _residual_`语句 2. 改用`dist=normal`做探索性分析	在试验设计阶段预估最小样本量：使用`proc power`计算Binary检验所需栏数，确保各组`prop_diarrhea`在0.2-0.8区间
边界主导	`Event`频数表中，某组`Frequency`为0或N；`proportion`分布直方图在0或1处有尖峰	`proc freq data=binomial_data; tables Y*N / list; run;`查看Y=0或Y=N的频数	1. 对Binomial数据，添加`random _residual_ / subject=pen*week` 2. 对Beta数据，增大扰动值至0.001	在数据收集阶段设置“腹泻触发阈值”：若连续3天score=0，暂停该栏数据采集，避免产生大量0值
模型过载	`Iteration History`中目标函数值震荡；`Covariance Parameter Estimates`表中多个参数为0或负值	`proc glimmix`后加`ods output ConvergenceStatus=converge;`检查`Reason`列	1. 移除最不重要的随机效应（如先删`block`） 2. 将`type=un`改为`type=vc`（方差成分）	采用“自下而上”建模：先拟合无随机效应模型，再逐个添加随机项，每步验证收敛性

提示：当所有解法失效时，我的终极手段是改用广义估计方程（GEE）。用proc genmod替代proc glimmix，指定repeated subject=pen / type=ar(1)。GEE不依赖分布假设，对收敛性更宽容，虽牺牲部分效率，但能给出稳健的边际效应估计。

5.2 残差诊断误区：为什么“看起来正常”反而是错的？

新手常犯的致命错误是：看到残差图呈“钟形”，就认为模型完美。Dr. Jacobs一针见血指出：“Using a Binary distribution, you cannot expect homogenous or normal error.” —— 这是统计学基本原理：Binary数据的残差必然服从二项分布，其方差=μ(1-μ)，天然异方差。正确诊断方法如下：

/* Binary残差诊断：应检查“残差 vs 预测概率”图 */ proc glimmix data=binary_trimmed plots(only)=effect(residuals); class treatment block pen; model event(event='1') = treatment / dist=binary link=logit solution; random block; random _residual_ / subject=pen type=ar(1); output out=resid_binary residual=r; lsmeans treatment / ilink cl; run; /* 解读要点： */ /* - 残差应在0附近随机散布，无趋势 */ /* - 靠近预测概率0或1的点，残差绝对值应较小（因方差小） */ /* - 若出现“漏斗形”（预测概率0.5处残差大，两端小），说明模型拟合良好 */

Beta残差的正确打开方式：
Beta分布的残差应近似正态，但需用Pearson残差而非原始残差。GLIMMIX中，plots(only)=effect(residuals)默认绘制Pearson残差。若Pearson残差图呈明显“U形”（两端高，中间低），说明beta分布的α/β参数未能充分捕捉比例分布的偏态，此时应尝试link=log或link=probit替代link=logit。

5.3 结果矛盾怎么办？——当Binary、Binomial、Beta给出不同结论时

这是高级用户才会遇到的困境。我经历过一个真实案例：Binary分析显示Treatment A显著优于B（p=0.02），Binomial分析无差异（p=0.15），Beta分析却显示A更差（p=0.04）。深入排查发现：

Binary的显著性源于少数高发栏：A组有2栏猪腹泻率90%，B组最高仅60%，但A组其余栏全为0%。Binary将“高发栏”权重过大，掩盖了整体趋势；
Binomial平滑了极端值：因按周聚合，A组高发栏的90%被稀释为每周平均，削弱了差异；
Beta捕捉了分布形态：A组比例分布双峰（0%和90%），B组单峰（40%），Beta的α/β参数检测到A组更高的离散度，故判定其“不稳定”。

决策树：

若临床关注事件发生风险（如“是否需要紧急用药”），以Binary结果为准；
若关注疾病负担总量（如“总治疗成本”），以Binomial结果为准；
若关注疾病动态稳定性（如“干预是否使病情更可控”），以Beta结果为准。
终极原则：没有“正确”的统计结果，只有“匹配科学问题”的结果。报告时必须明确声明：“本分析采用Binary分布，旨在评估Treatment对腹泻发生率的影响”，而非模糊地说“分析显示Treatment有效”。

6. 从实验室到论文：结果呈现与可视化的一线经验

6.1 图表制作的黄金三原则：信息密度、临床可读性、统计严谨性

在向兽医同行汇报时，我坚持三个图表原则：

原则一：拒绝“学术黑箱”图——所有图表必须能让不懂SAS的人看懂横纵轴含义。例如，Binary的lsmeans图，Y轴必须标注“Estimated Probability of Diarrhea (%)”，而非“Logit Scale Estimate”；
原则二：误差线必须匹配分布——Binary/Binomial用Wald置信区间（cl选项输出），Beta用profile likelihood置信区间（需nloptions tech=nms启用）；
原则三：时间趋势图必须标注聚合粒度——若用周数据，图中必须注明“Weekly Aggregated Data (N=7 days/week)”，避免读者误以为是日数据。

实操模板（SAS GTL代码）：

/* Binary结果图：各处理组腹泻发生率 */ proc template; define statgraph binary_plot; begingraph; entrytitle "Diarrhea Incidence by Treatment"; layout overlay / yaxisopts=(label="Probability of Diarrhea" linearopts=(viewmin=0 viewmax=1)); scatterplot x=treatment y=estimate / yerrorlower=lcl ilink yerrorupper=ucl ilink name="scatter" markerattrs=(size=10); discretelegend "scatter" / title="Treatment"; endlayout; endgraph; end; run; proc sgrender data=lsmeans_binary template=binary_plot; run;

6.2 论文写作的隐藏陷阱：如何描述方法而不暴露统计短板

审稿人最常质疑的是方法学描述。我总结了安全表述的三句话模板：

描述Binary：“We modeled diarrhea occurrence (score ≥2) as a binary outcome using generalized linear mixed models (GLIMMIX procedure, SAS v9.4), with treatment as fixed effect and batch (block) and pen-specific autocorrelation as random effects.”
描述Binomial：“Weekly diarrhea incidence was calculated as the proportion of diarrheic days per pen, modeled using a binomial distribution with logit link to account for the bounded nature of proportion data.”
描述Beta：“To characterize the continuous severity of diarrhea, we modeled the weekly proportion of diarrheic days using a beta distribution, which explicitly accommodates the heteroscedasticity inherent in bounded proportion data.”

绝对避免的表述：

❌ “We chose Beta distribution because it gave better fit”（暴露模型比较，易被质疑p-hacking）；
✅ “Given our objective to quantify diarrheal severity as a continuous measure, the beta distribution was selected a priori for its theoretical appropriateness in modeling bounded proportions.”（强调先验合理性）。