蛋鸡死亡发生率建模：从零膨胀计数到产业级预警

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_intervalDate记录日期日粒度所有时间序列分析的锚点必须连续无缺失。需检查是否为自然日非生产日历lcn_idFactor鸡舍IDLocation ID捕获鸡舍间固有差异建筑、设备、管理风格必须为分类变量避免编码为数字如1436导致模型误判为连续变量animals_startNumeric当日晨间存栏量头计算死亡率的分母是“风险人口”的基线是核心分母必须严格校验。常见错误将“转入数”误加或“转出数”未扣减animals_deadNumeric当日死亡数头模型的因变量Response Variable是核心分子必须为整数。需重点筛查异常值如单日死亡100只需人工复核t_min_inNumeric当日舍内最低温度℃冷应激指标尤其在夜间或通风过度时与t_max_in共同构成温差是热应激的关键代理变量t_max_inNumeric当日舍内最高温度℃热应激核心指标与死亡率呈强正相关需与t_min_in配对使用单独看意义有限t_act_outNumeric当日舍外实际温度℃衡量环控系统负荷反映外部环境压力是评估环控效能的外部参照其他如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的分布直方图bins20足够并计算其零比例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_in0) %% filter(t_max_in0)这看似是过滤负值实则是规避了更深层的问题负温在鸡舍中本就不该出现出现即意味着传感器严重漂移或数据传输错误此时应追溯原始传感器日志而非简单粗暴地删除。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 ~ ..., familybinomial)并将event定义为ifelse(animals_dead0, 1, 0)。这是一个极具迷惑性的起点因为它看起来无比自然有死亡就是1没死亡就是0不就是经典的“是/否”分类问题吗然而这恰恰是初学者最容易踏入的统计学陷阱。让我用一个生活化的例子来解释假设你是一家保险公司的精算师要为一辆车定价。你有两个选择A. 统计这辆车在过去一年里“是否出过事故”是/否B. 统计这辆车在过去一年里“出了几次事故”0次、1次、2次…。显然B选项包含了A选项的全部信息如果出了2次那肯定“是”出过还额外提供了事故频次的强度信息。二项模型只用了A选项它把“死亡1只”和“死亡50只”完全等同起来都记为event1。这在统计上叫做信息的巨大损失。更致命的是二项模型的数学假设与蛋鸡死亡的现实严重冲突。二项分布要求每一次“试验”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(Yk) (λ^k * e^{-λ}) / k!简洁而强大。在R中一行代码即可拟合glm(animals_dead ~ temp_range t_max_in offset(log(animals_start)), familypoisson)。这里的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日窗口则是参考了鸡只生理恢复的典型周期。这些都不是凭空而来而是十年一线经验的结晶。5.3 模型拟合、