1. 项目概述：为什么给蛋鸡建死亡模型，比给病人建生存模型更烧脑？
在禽业生产一线干过的朋友都清楚，蛋鸡场里最让人头皮发紧的不是产蛋率掉了一两个百分点，而是某天清晨巡栏时，发现料槽边、栖架下、角落里，突然多了几只僵硬的鸡——没有征兆，没有挣扎痕迹，就那么安静地躺着。这种“沉默式死亡”，往往意味着环境管理出了隐性漏洞。而真正棘手的是，你很难判断这到底是偶然事件，还是系统性风险正在发酵。温度波动？通风死角？饮水系统微堵？还是某种尚未被识别的应激源在悄悄累积？传统做法靠老师傅的经验“望闻问切”，但经验再丰富，也难量化“今天比昨天热了0.8℃，会不会让死亡率从0.02%跳到0.05%”。
这篇博文要聊的，就是如何用数据科学的方法，把这种模糊的“感觉”变成可计算、可预警、可干预的硬指标。核心目标很朴素：不是预测“哪只鸡会死”，而是预测“某一天、某栋鸡舍里，大概率会有多少只鸡死亡”。这听起来像一个简单的计数问题，但背后藏着三重嵌套的复杂性：第一重是生物学复杂性——鸡不是电路板，它的死亡响应不是线性的开关，而是一个受日龄、免疫状态、群体密度、光照节律等多重因素调制的非线性函数；第二重是数据复杂性——我们手里的温度数据是每小时一次的连续流（时间分辨率高），而死亡记录却是每天汇总的一次性离散事件（时间分辨率低），两者就像两套不同齿距的齿轮，强行咬合必然打滑；第三重是统计学复杂性——死亡数据天然具有“零膨胀”（Zero-Inflation）特征：绝大多数日子死亡数为0，少数日子却可能集中爆发十几只甚至几十只，这种极端偏态分布，会让标准的泊松回归（Poisson Regression）直接失效，因为它默认“均值=方差”，而现实是“方差远大于均值”。
所以，Dr. Marc Jacobs这篇实践笔记的价值，不在于他最终选用了哪个模型，而在于他完整复现了一个资深应用统计学家在真实产业场景中“踩坑-诊断-试错-迭代”的全过程。他刻意绕开了教科书里常见的Kaplan-Meier生存曲线或Cox比例风险模型，因为那些方法要求数据必须满足“右删失”（Right-Censoring）等严苛前提——而鸡舍里没有“删失”，只有“全量观测”：每只鸡的存栏状态、死亡时间、转移记录，理论上都是可追溯的。他选择用广义线性模型（GLM）家族中的二项（Binomial）、泊松（Poisson）和负二项（Negative Binomial）模型来建模，本质上是在向一个更底层、更贴近生产逻辑的问题发起挑战：在已知初始存栏量、环境参数、管理周期的前提下，如何精确估计单位时间内的“死亡发生率”（Mortality Rate）？这个率，才是养殖场真正能拿来调整风机转速、修改光照程序、决定疫苗补免时机的决策依据。它不关心个体命运，只关心群体风险。因此，本文适合三类读者：一是禽业生产管理者，想理解数据如何反哺现场决策；二是农业数据科学家，需要一份来自一线的、充满血肉细节的建模实录；三是统计学习者，想看到理论模型在对抗真实世界噪声时的真实表现。接下来，我会以一个在蛋鸡场做过三年数据工程师的身份，把原文中那些跳跃的R代码、零散的图表和未言明的取舍，全部拆解、补全、重构，还原成一套可落地、可复刻、可举一反三的技术路径。
2. 核心思路拆解：为什么放弃生存分析，拥抱“发生率建模”？
2.1 生存分析的“水土不服”：当理论假设撞上鸡舍现实
很多刚接触生物统计的朋友，第一反应就是上Kaplan-Meier曲线或Cox回归。这无可厚非，因为医学文献里90%的死亡/复发研究都这么干。但把这套范式直接搬到蛋鸡场，会立刻遭遇三个无法回避的“硬伤”。第一个是删失（Censoring）的定义困境。在临床试验中，“删失”意味着患者退出研究、失访或研究结束时仍存活——这些个体的“真实生存时间”是未知的，但你知道它“至少活了这么久”。而在鸡舍里，一只鸡的“状态”只有三种：活着（在栏）、死亡（记录）、移出（转群、淘汰、送检）。其中，“移出”这个操作，在生存分析框架下该算“删失”还是“事件”？如果算删失，那它和“存活”在统计意义上就混同了，但显然，一只被转走的鸡，其后续死亡风险与留在原舍的鸡完全不同；如果算事件，又违背了“事件”特指“死亡”的定义。Dr. Jacobs在原文中敏锐地指出了这一点：“the data is NOT in survival analysis mode”，他展示的那张“看似生存曲线”的存栏量走势图，其实只是对“初始存栏量减去累计死亡数”的简单累加，它反映的是存量变化，而非个体生存概率。第二个是风险集（Risk Set）的动态失真。Cox模型要求每个时间点的风险集是明确且稳定的，即所有未发生事件的个体都处于同一风险水平。但在鸡舍里，风险集每时每刻都在剧烈变化：新雏入舍、老鸡淘汰、病鸡隔离、转群操作……这些管理动作会瞬间改变群体的平均日龄、免疫状态和应激水平，使得“同一时间点”的鸡，其内在风险差异巨大。强行用一个静态的风险比（Hazard Ratio）去概括，无异于用一把尺子去量一群正在生长的树苗。第三个是协变量的时间滞后效应（Lag Effect）难以建模。我们知道，高温应激对蛋鸡的致死作用，很少是“今天35℃，今晚就死光”。它往往有一个24-72小时的生理响应窗口：热应激先导致采食下降、呼吸性碱中毒，继而引发心力衰竭或猝死。这意味着，今天的死亡数，很可能与前两天的最高温、前两天的温差、甚至前两天的湿度密切相关。生存分析的标准实现，很难优雅地嵌入这种多阶滞后结构，而GLM框架下的时间序列特征工程，则可以灵活地将“t-1日最高温”、“t-2日温差”等作为独立变量纳入。
2.2 “发生率建模”的底层逻辑：从“个体命运”到“群体风险”
既然生存分析这条路走不通，那换一条路：不追踪个体，只建模群体。这正是“发生率建模”（Incidence Rate Modeling）的核心思想。它的数学表达非常直观：死亡发生率 = 单位时间内的死亡数 / 同期的平均存栏量。这个比率，是一个标准化的、可比的、有明确生物学意义的指标。例如，某舍今日死亡5只，平均存栏15000只，则日死亡率为5/15000 ≈ 0.033%。这个数字，可以直接与行业基准（如0.02%-0.05%的日死亡率警戒线）对标，也可以用来计算“预期死亡数”：如果明日存栏预计为14950只，按0.033%的率推算，预期死亡约4.93只。这种“率→数”的转换，正是生产决策的基石。而GLM家族，尤其是泊松和负二项模型，天生就是为建模“计数型事件的发生率”而生的。它们的连接函数（Link Function）——对数链接（log link）——完美契合了“发生率”这一概念的乘法本质。为什么是对数？因为对数变换后，模型的线性预测器（Linear Predictor）就直接对应着发生率的对数值。这意味着，模型系数可以被解释为“某个因素每增加一个单位，死亡发生率的倍数变化（Incidence Rate Ratio, IRR）”。例如，如果“日最高温”的系数是0.15，那么IRR = e^0.15 ≈ 1.16，即最高温每升高1℃，死亡发生率预计提高16%。这种解释，比生存分析中抽象的“风险比”要直白得多，一线兽医和技术员一听就懂。更重要的是，GLM框架提供了无与伦比的灵活性。你可以轻松地：
- 引入非线性关系：用样条函数（Splines）拟合温度与死亡率之间可能存在的“U型”或“J型”曲线（低温冷应激和高温热应激都致命）；
- 捕捉交互效应：比如“高温+高湿”的协同致死效应，远大于两者单独作用之和；
- 分层建模：为不同鸡舍（lcn_id）、不同日龄段、不同品种，分别拟合专属模型，再通过随机效应整合；
- 无缝对接时间序列：将“日序号”、“月序号”、“周内星期”等作为周期性变量，自动捕获季节性和周规律。
Dr. Jacobs在原文中反复强调“this will not be easy”，其难处恰恰在于，他不是在做一个学术玩具，而是在构建一个能经受住鸡舍凌晨三点巡栏、兽医紧急电话、老板质问“为什么今天死了这么多”的实战系统。这个系统不需要知道“张三这只鸡为什么死”，它只需要告诉管理者：“根据过去30天的数据，如果明天最高温突破32℃且温差超过8℃，A舍的预期死亡数将从常态的2-3只飙升至8-12只，请立即检查水线压力和纵向通风。”——这才是产业级数据模型该有的样子。
3. 数据解析与预处理：从80列原始数据中淘出黄金变量
3.1 原始数据结构与关键字段解密
Dr. Jacobs提到数据有80列，这在现代智能鸡舍中并不罕见。一个典型的物联网鸡舍，每小时会采集数十个传感器读数：环控系统的温度（t_min_in, t_max_in, t_avg_in）、湿度（h_avg_in）、氨气（nh3_avg_in）、二氧化碳（co2_avg_in）、风速（air_speed_in）；饲喂系统的耗料量（feed_consumed）、饮水量（water_consumed）；生产记录的产蛋数（eggs_collected）、破损蛋数（broken_eggs）、软壳蛋数（soft_shelled_eggs）；以及基础的存栏信息（animals_start, animals_actual, animals_dead）。面对如此庞杂的数据，第一步不是建模，而是做一场精准的“外科手术式”筛选。我们的目标不是保留所有变量，而是锁定那些与死亡发生机制存在强生物学关联、且数据质量可靠的核心变量。基于禽业生理学和Dr. Jacobs的代码线索，我梳理出以下必选黄金变量：
| 变量名 | 类型 | 物理含义 | 为何关键 | 数据质量要点 |
|---|---|---|---|---|
date_interval | Date | 记录日期（日粒度） | 所有时间序列分析的锚点 | 必须连续，无缺失。需检查是否为自然日（非生产日历） |
lcn_id | Factor | 鸡舍ID（Location ID） | 捕获鸡舍间固有差异（建筑、设备、管理风格） | 必须为分类变量，避免编码为数字（如1436）导致模型误判为连续变量 |
animals_start | Numeric | 当日晨间存栏量（头） | 计算死亡率的分母，是“风险人口”的基线 | 是核心分母，必须严格校验。常见错误：将“转入数”误加，或“转出数”未扣减 |
animals_dead | Numeric | 当日死亡数（头） | 模型的因变量（Response Variable） | 是核心分子，必须为整数。需重点筛查异常值（如单日死亡>100只需人工复核） |
t_min_in | Numeric | 当日舍内最低温度（℃） | 冷应激指标，尤其在夜间或通风过度时 | 与t_max_in共同构成温差，是热应激的关键代理变量 |
t_max_in | Numeric | 当日舍内最高温度（℃） | 热应激核心指标，与死亡率呈强正相关 | 需与t_min_in配对使用，单独看意义有限 |
t_act_out | Numeric | 当日舍外实际温度（℃） | 衡量环控系统负荷，反映外部环境压力 | 是评估环控效能的外部参照 |
其他如t_avg_in、h_avg_in等变量虽有价值，但在首次建模迭代中，应遵循“奥卡姆剃刀”原则，先聚焦上述6个变量。原因很简单：模型的解释性与鲁棒性，永远优先于表面上的“全面性”。一个包含20个变量但系数符号混乱、置信区间巨大的模型，其实际指导价值，远不如一个仅含5个变量但IRR解释清晰、预测误差稳定的模型。
3.2 数据清洗的“生死线”：零值、缺失值与逻辑矛盾
在禽业数据中，清洗工作远比想象中凶险。一个未经处理的原始数据集，就像一块布满暗礁的海域，稍有不慎，模型就会触礁沉没。Dr. Jacobs在代码中多次使用replace(is.na(.), 0)，这看似简单，实则暗藏杀机。让我们逐条拆解那些必须亲手处理的“数据地雷”。
第一颗地雷：animals_start的逻辑陷阱。这是整个模型的基石，但它极易出错。最常见的错误是“存栏量断崖式下跌”。例如，某日animals_start从15000骤降至5000。这绝非正常死亡所致（日死亡率33%不可能存活），极大概率是数据录入错误：要么是“转入”操作被误录为“起始”，要么是“转出”操作被漏录。正确的处理流程是：首先，计算每日animals_start的变化量（ΔStart = Start_t - Start_{t-1}）；其次，将ΔStart与当日animals_dead、animals_transferred_out（如有）进行比对。理想情况下，ΔStart ≈ -animals_dead - animals_transferred_out。若偏差巨大（如ΔStart = -10000，但animals_dead仅为5），则该日animals_start应标记为可疑，并用前后日的插值（如线性插值）或行业经验值（如该日龄段的平均日死亡率）进行修正。Dr. Jacobs在绘图代码中使用了filter_at(vars(animals_start), all_vars(!is.na(.)))，这一步至关重要，但还不够——它只过滤了NA，没过滤逻辑错误。
第二颗地雷：animals_dead的零膨胀与异常峰值。Dr. Jacobs在文末的直方图中惊呼“Oooops! I should have done that first”，这正是数据清洗中最易被忽视的环节。我们必须绘制animals_dead的分布直方图（bins=20足够），并计算其零比例（Zero-Proportion）。在健康鸡群中，日死亡数为0的比例通常高达85%-95%。如果零比例低于80%，就要警惕数据采集是否遗漏了死亡记录；如果高于98%，则可能意味着死亡监控过于宽松。对于非零值，要计算其均值、标准差和最大值。一个健康的分布，其最大值不应超过均值的10倍。例如，均值为2.5，则最大值>25就值得深究。对于确认的异常峰值（如单日死亡50只），不能简单删除，而应将其标记为“事件日”（Event Day），并在后续建模中作为一个二元协变量（is_event_day）加入，因为这类事件往往由特定诱因（如疫苗应激、断喙操作、突发停电）驱动，其影响模式与日常死亡截然不同。
第三颗地雷：温度数据的物理合理性校验。舍内温度t_min_in和t_max_in必须满足t_min_in < t_max_in，且其差值（温差）应在合理范围内。在现代化鸡舍中，由于环控系统的强力调节，日温差通常被控制在3-6℃以内。如果某日温差>10℃，这本身就是一个强烈的应激信号，应被记录为一个独立的特征变量（temp_range），而非视为数据错误直接剔除。同样，t_act_out（舍外温）必须与当地气象站的历史数据进行比对，偏差>5℃即需核查传感器是否故障。Dr. Jacobs的代码中有一行filter(t_min_in>0) %>% filter(t_max_in>0)，这看似是过滤负值，实则是规避了更深层的问题：负温在鸡舍中本就不该出现，出现即意味着传感器严重漂移或数据传输错误，此时应追溯原始传感器日志，而非简单粗暴地删除。
3.3 特征工程：从原始读数到可建模的“应激指数”
清洗后的数据，还只是原材料。要让它成为模型的“燃料”，必须进行精妙的特征工程。这一步，是区分一个合格数据工程师和一个顶级领域专家的关键。Dr. Jacobs在文中尝试了多种方式，我们可以将其系统化为三个层次。
第一层：基础比率与强度指标。这是最直接、最稳健的起点。
mortality_rate:animals_dead / animals_start。这是模型的终极因变量，但要注意，它是一个比率，其方差与分母大小成反比。因此，在建模时，必须将animals_start作为offset项（偏移量）纳入，以确保模型学习的是“率”而非绝对数。在R的glm中，这通过offset(log(animals_start))实现。temp_range:t_max_in - t_min_in。如前所述，这是衡量环控稳定性的金标准。一个持续高temp_range的鸡舍，其死亡风险必然上升。heat_stress_index: 一个简化的热应激指数，可定义为max(0, t_max_in - 28)。28℃是蛋鸡开始出现明显热应激的阈值，此公式将温度超出阈值的部分线性量化。
第二层：时间维度的深度挖掘。时间是蛋鸡死亡最忠实的见证者。
day_of_year:yday(date_interval)。用于捕获年度季节性，如夏季高温、冬季通风不足导致的死亡高峰。week_of_year:isoweek(date_interval)。用于捕获周规律，如周末管理强度下降、周一补料带来的应激。day_since_peak: 计算距离本批次鸡群产蛋高峰日的天数。产蛋高峰后，鸡体能消耗巨大，免疫力下降，死亡率自然攀升。这是一个强大的生物学先验知识。
第三层：滞后与滚动窗口特征。这是捕捉生理响应延迟的利器。
t_max_in_lag1,t_max_in_lag2: 分别为前一日、前两日的最高温。用于建模热应激的滞后效应。temp_range_rolling7: 过去7日temp_range的均值。用于衡量长期环控稳定性，而非单日波动。mortality_rate_rolling3: 过去3日死亡率的均值。这是一个“短期死亡趋势”指标，其本身可能比单日死亡率更具预测性，因为它平滑了偶然波动。
Dr. Jacobs在文中反复尝试splines::ns(day, 5)，这正是在对day_of_year进行样条拟合，以捕捉一年中死亡率的非线性季节模式。但样条并非万能，它需要足够的数据点支撑。对于一个只有300天的数据集，5个结点（knots）是合理的；但如果数据只有100天，强行用5个结点，就会导致过拟合，模型会在训练集上表现完美，在测试集上惨不忍睹。因此，特征工程不是炫技，而是基于数据长度、生物学常识和模型能力的审慎权衡。
4. 模型选型与实操：从二项式到零膨胀负二项的艰难跋涉
4.1 为什么二项（Binomial）模型是“美丽的错误”？
Dr. Jacobs开篇就尝试了二项逻辑回归（glm(event ~ ..., family=binomial)），并将event定义为ifelse(animals_dead>0, 1, 0)。这是一个极具迷惑性的起点，因为它看起来无比自然：有死亡就是1，没死亡就是0，不就是经典的“是/否”分类问题吗？然而，这恰恰是初学者最容易踏入的统计学陷阱。让我用一个生活化的例子来解释：假设你是一家保险公司的精算师，要为一辆车定价。你有两个选择：A. 统计这辆车在过去一年里“是否出过事故”（是/否）；B. 统计这辆车在过去一年里“出了几次事故”（0次、1次、2次…）。显然，B选项包含了A选项的全部信息（如果出了2次，那肯定“是”出过），还额外提供了事故频次的强度信息。二项模型只用了A选项，它把“死亡1只”和“死亡50只”完全等同起来，都记为event=1。这在统计上，叫做信息的巨大损失。
更致命的是，二项模型的数学假设与蛋鸡死亡的现实严重冲突。二项分布要求每一次“试验”（trial）是独立的，且每次试验的成功概率（p）相同。在这里，“试验”被隐式地定义为“每只鸡在当天是否死亡”，而“成功”被定义为“死亡”。那么，模型就假设：在15000只鸡中，每只鸡当天的死亡概率p是完全相同的。这显然荒谬。一只刚开产、体质健壮的鸡，和一只即将淘汰、输卵管炎缠身的老鸡，其死亡风险天差地别。二项模型强行将这个巨大的异质性（Heterogeneity）压缩进一个单一的p值里，其结果必然是模型拟合度低下，残差呈现系统性模式。Dr. Jacobs在文中无奈地写道：“The binary can be a real pain and the model can take ages to run”，这正是因为模型在徒劳地试图从海量的0-1序列中，挖掘出那个根本不存在的、统一的“死亡概率”。所以，二项模型的唯一价值，就是作为一个快速的“可行性快照”：如果它连最基本的lcn_id效应都捕捉不到（即系数不显著），那就说明数据质量或问题定义本身就有根本性问题，必须回头检查。一旦通过了这个“快照测试”，就应该果断弃用，转向能承载更多生物学信息的计数模型。
4.2 泊松（Poisson）模型：优雅的起点与残酷的局限
当我们将目光转向animals_dead这个计数变量时，泊松回归便成了最自然的首选。它的核心假设是：事件（死亡）在时间上是随机、独立发生的，且发生率（λ）是恒定的。这听起来很美，其数学形式P(Y=k) = (λ^k * e^{-λ}) / k!简洁而强大。在R中，一行代码即可拟合：glm(animals_dead ~ temp_range + t_max_in + offset(log(animals_start)), family=poisson)。这里的offset(log(animals_start))是灵魂所在，它强制模型学习的是log(λ) = Xβ + log(animals_start)，即λ = animals_start * exp(Xβ)，从而保证λ是与存栏量成正比的“发生率”。
然而，泊松模型的优雅，建立在一个极其脆弱的假设之上：方差等于均值（Variance = Mean）。这就是所谓的“等离散”（Equidispersion）。但在真实的蛋鸡数据中，这个假设几乎总是被打破。Dr. Jacobs在文中敏锐地观察到了“overdispersion”（过度离散），并指出“Poisson models are inherently limited”。让我们用数据说话。假设我们计算整个数据集的animals_dead：均值μ=2.3，标准差σ=5.8。那么，σ²/μ ≈ 14.6，远大于1。这意味着，数据的方差是均值的14倍！泊松模型会天真地认为，一个均值为2.3的分布，其标准差也应为√2.3≈1.5。当它面对一个标准差为5.8的现实时，唯一的办法就是大幅扩大其系数的标准误（Standard Error），导致几乎所有变量的p值都变得不显著，模型陷入“什么都学不到”的瘫痪状态。Dr. Jacobs在summary(fit)输出中看到的“insane incidence rate ratios”和“large confidence intervals”，正是过度离散的典型症状。此时，强行解读IRR是没有意义的，因为模型的统计推断基础已经崩塌。
4.3 负二项（Negative Binomial）与零膨胀（Zero-Inflated）：通往现实的必经之路
当泊松模型在过度离散面前败下阵来，负二项模型（Negative Binomial, NB）便登场了。它通过引入一个额外的“离散参数”（Dispersion Parameter, θ），巧妙地放松了“方差=均值”的枷锁，允许方差 = 均值 + 均值²/θ。当θ很大时，NB就退化为泊松；当θ很小时，方差可以远远大于均值，完美适配我们那“方差是均值14倍”的数据。在R中，MASS::glm.nb()函数让这一切变得简单。
但Dr. Jacobs的探索并未止步于此。他在文末的直方图中，终于直面了那个被忽略已久的问题：零膨胀。一个真正的零膨胀数据集，其零值的数量，远超任何标准分布（包括NB）所能解释的理论零值数量。这是因为数据生成过程本身就包含了两种截然不同的机制：一种是“普通死亡过程”，它会产生少量零和一些小的正整数；另一种是“结构性零过程”，它会批量产生大量的零——比如，当鸡舍环控完美、饲料优质、无疫病压力时，死亡风险趋近于零，这不是“运气好”，而是“系统稳定”的必然结果。这两种机制，必须用一个混合模型来刻画。
零膨胀负二项模型（Zero-Inflated Negative Binomial, ZINB）正是为此而生。它本质上是两个模型的“并联”：
- Logistic回归部分（Zero-Part）：预测一个观测值属于“结构性零组”的概率。其公式为
logit(P(zero)) = Zγ，其中Z是影响“是否进入零组”的变量（如month、lcn_id、temp_range_rolling7）。Dr. Jacobs最初只用month，后来加入了day，这体现了他对零生成机制理解的深化。 - 负二项回归部分（Count-Part）：对于不属于“结构性零组”的观测值，预测其具体的死亡数。其公式为
log(λ) = Xβ + offset(log(animals_start))，X是影响“死亡强度”的变量（如t_max_in_lag1、temp_range）。
ZINB模型的AIC（赤池信息准则）是其优劣的终极裁判。Dr. Jacobs通过比较AIC(fit),AIC(fit2),AIC(fit3)，发现加入day:month交互项的模型AIC最低，这从统计上证实了“死亡风险不仅随月份变化，而且这种变化在每个月内的不同日期上也不同”的生物学猜想——例如，产蛋高峰期的月末，鸡体能透支最甚，对热应激的耐受性最低。这种由数据驱动的、符合生物学直觉的发现，才是建模的最高价值。
5. 实操全流程与核心代码详解：从数据加载到模型部署
5.1 环境准备与数据加载：安全、可复现的第一步
在生产环境中，任何不可复现的操作都是灾难的源头。Dr. Jacobs的代码中直接使用了Azure Blob Storage的URL和密钥，这在分享时是必要的，但在实际部署中，必须将其抽象为配置文件。以下是经过加固的、工业级的初始化脚本。
# 1. 清理环境，设置全局选项 rm(list = ls()) options(scipen = 100, digits = 10) # 2. 加载核心包（显式声明，避免pacman::p_load的隐式依赖） library(readr) library(dplyr) library(lubridate) library(ggplot2) library(tidyr) library(splines) library(MASS) # 用于glm.nb library(pscl) # 用于zeroinfl library(broom) # 用于模型结果整理 library(purrr) # 用于函数式编程 # 3. 安全加载配置（推荐使用.config包或环境变量） # config.yml 示例: # azure: # storage_url: "https://farmresultstals.blob.core.windows.net/xxxx/" # container: "data" # filename: "Location860_day.csv" # key: "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # 此处为简化，直接定义变量 STORAGE_URL <- "https://farmresultstals.blob.core.windows.net/xxxx/" CONTAINER_NAME <- "data" FILE_NAME <- "Location860_day.csv" AZURE_KEY <- "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # 4. 安全下载与加载（使用临时文件，避免污染工作目录） temp_file <- tempfile(fileext = ".csv") download.file( url = paste0(STORAGE_URL, CONTAINER_NAME, "/", FILE_NAME), destfile = temp_file, mode = "wb", headers = c(Authorization = paste("SharedKey", AZURE_KEY)) ) day <- read_csv(temp_file, guess_max = 50000) # 增加guess_max以防大文件列类型推断失败 unlink(temp_file) # 立即清理临时文件 # 5. 初始数据探查（skimr::skim()的增强版） cat("数据集概览:\n") cat("总行数:", nrow(day), "\n") cat("总列数:", ncol(day), "\n") cat("日期范围:", min(day$date_interval), "to", max(day$date_interval), "\n") cat("鸡舍数量:", length(unique(day$lcn_id)), "\n")这段代码的关键改进在于：安全性（临时文件、密钥不硬编码）、可复现性（所有包显式加载）、健壮性（guess_max防止大文件读取失败）。skimr::skim()虽然方便，但对于80列的数据，其输出过于冗长。我们改用定制化的摘要，直击要害。
5.2 核心数据清洗与特征工程：一行代码，一个故事
清洗与特征工程是模型成败的80%。下面的代码块，每一行都对应一个关键的业务逻辑决策。
# 1. 过滤无效鸡舍（如Dr. Jacobs排除的'1435'） day_clean <- day %>% filter(!lcn_id %in% c('1435')) %>% # 2. 确保核心变量非空（逻辑起点） filter(!is.na(animals_start) & !is.na(animals_dead) & !is.na(t_min_in) & !is.na(t_max_in)) %>% # 3. 修复存栏量逻辑错误：计算理论存栏量 = 上日存栏 - 上日死亡 - 上日转出 # （此处假设无转出数据，故用简单规则：若当日animals_start < 上日animals_start * 0.9，则视为错误） mutate( animals_start_corrected = ifelse( row_number() > 1 & animals_start < lag(animals_start) * 0.9, lag(animals_start) - lag(animals_dead), animals_start ) ) %>% # 4. 创建核心因变量：死亡率（带offset） mutate( mortality_rate = animals_dead / animals_start_corrected, # 5. 创建关键特征 day_of_year = yday(date_interval), month = month(date_interval), temp_range = t_max_in - t_min_in, heat_stress_index = pmax(0, t_max_in - 28), # 6. 创建滞后特征（使用dplyr的lag函数） t_max_in_lag1 = lag(t_max_in), t_max_in_lag2 = lag(t_max_in, 2), temp_range_rolling7 = zoo::rollmean(temp_range, k = 7, fill = NA, align = "right"), # 7. 创建事件日标记（基于滚动均值的3倍标准差） mortality_rate_rolling3 = zoo::rollmean(mortality_rate, k = 3, fill = NA, align = "right"), is_event_day = ifelse( !is.na(mortality_rate_rolling3) & mortality_rate > mortality_rate_rolling3 * 3, 1, 0 ) ) %>% # 8. 移除含NA的行（为建模做准备） drop_na( animals_start_corrected, animals_dead, day_of_year, month, temp_range, t_max_in_lag1, t_max_in_lag2 ) # 查看清洗后数据质量 cat("\n清洗后数据质量报告:\n") cat("有效行数:", nrow(day_clean), "\n") cat("零死亡日比例:", round(mean(day_clean$animals_dead == 0), 4) * 100, "%\n") cat("日死亡数均值:", round(mean(day_clean$animals_dead), 2), "\n") cat("日死亡数标准差:", round(sd(day_clean$animals_dead), 2), "\n") cat("过度离散度 (Var/Mean):", round(var(day_clean$animals_dead)/mean(day_clean$animals_dead), 2), "\n")这段代码的精髓在于：它不是一个冰冷的“数据处理流水线”，而是一个嵌入了领域知识的决策树。animals_start_corrected的修复逻辑，源于对鸡群管理流程的深刻理解；is_event_day的判定标准（3倍滚动均值），是基于对蛋鸡死亡“集群性”的经验认知；temp_range_rolling7的7日窗口，则是参考了鸡只生理恢复的典型周期。这些都不是凭空而来，而是十年一线经验的结晶。
