之前有个项目要处理一份个人消费问卷样本量三千多因变量是“过去三个月网购花费”。数据里四成左右的人填了0剩下的人金额从十几块到几万块不等。组里一开始直接跑了OLS系数看着显著但预测值出现大量负数领导看了一眼就否了。后来换成tobit模型才把问题说清楚。今天就专门聊一个组合VGAM包 tobit模型把审查回归的原理、R实现、边际效应和实际项目里踩过的坑一次讲透。如果你要处理工资、支出、理赔额这类“下限或上限处堆积”的数据这篇应该能帮上忙。1. tobit模型的核心机制0堆不是你想删掉就能删掉1.1 一个典型的“零堆积”数据长什么样受限因变量在真实数据里太常见了。消费者支出问卷里没买过的人填0劳动收入数据里失业的人收入是0保险理赔记录里没出险的人理赔额是0。这些0和“真正的0”不一样它们背后对应的是一个本应为负但被观测机制掩盖的潜在值。假设一个人的真实消费意愿是 y*可以是负数但问卷只能记录真实的消费金额不能记录“负消费”。于是我们看到的 y 是这样生成的当 y* ≤ 0 时y 0当 y* 0 时y y*。这就是左侧审查left censoring。阈值0并不是随机出现的而是数据收集机制造成的。处理这种数据第一反应通常是两种要么把0当成普通数值直接跑OLS要么嫌0碍事删掉只对正数部分回归。这两种做法都会出问题。1.2 潜变量与审查机制的似然表达tobit模型把观测结果拆成两部分来处理。假设y* Xβ εε ~ N(0, σ²)观测到的 y 满足当 y* ≤ L 时y L当 y* L 时y y*。L 就是审查点censoring point比如0。对第 i 个观测似然贡献要分情况写。如果 y_i 正好落在审查点 L 上我们能确定的是 y*_i ≤ L所以贡献是累积概率P(y* ≤ L) Φ((L - Xβ) / σ)如果 y_i L我们观测到了具体数值贡献就是正态密度f(y_i) (1/σ) · φ((y_i - Xβ) / σ)整个模型的似然函数就是这两个部分的乘积。极大似然估计会同时估计 β 和 σ因为审查概率里同时出现了 β 和 σ不能像OLS那样把 σ 消掉。这也是tobit和普通线性回归在估计机制上的根本差异。1.3 为什么OLS和“删0再回归”都会翻车OLS把所有0当成真实的连续观测相当于假设 E(y|x) Xβ。但数据的真实期望是E(y|x) P(y*0|x) · E(y*|y*0, x)这是一个非线性函数用一条直线去拟合斜率必然被拉偏。更直观的问题是OLS完全无法保证预测值非负大量预测值掉到0以下业务上无法解释。把0删掉只对正数部分回归问题更隐蔽。y* 的条件分布原本是正态分布但条件在 y* 0 之后分布就变成了截断正态。截断正态的均值不再是 Xβ而是 Xβ σλ其中 λ 是逆米尔斯比率。直接用OLS拟合截断后的样本遗漏了这一项就会造成遗漏变量偏差系数同样不靠谱。所以tobit模型不是“换一个回归方法”的问题而是模型设定上要和数据的生成机制保持一致。2. 认识VGAM包里的tobit家族Lower/Upper和参数结构2.1 VGAM与AER::tobit的定位差异很多人第一次做tobit用的是AER::tobit()这个函数确实方便一行公式就能出结果。它适合快速交付但如果你的模型需求稍微复杂一点比如要做双审查、要在多个分布族之间切换对比或者希望在一个统一的框架下理解不同受限模型那么VGAM包是更好的选择。VGAM 的全称是 Vector Generalized Linear/Additive Models核心思想是把广义线性模型的多参数分布族统一起来。vglm()是它的主拟合函数family 参数可以是 tobit、probit、logit、负二项、零膨胀泊松等等。tobit 只是其中的一个族整个框架是通用的。简单说AER::tobit()是“专用小工具”VGAM::vglm()是“带抽屉的工具箱”。如果你要做模型比较、快速扩展VGAM 的灵活度明显更高。2.2 tobit()函数签名与审查边界设置VGAM 里拟合 tobit 的标准写法是library(VGAM) fit - vglm(y ~ x1 x2, tobit(Lower 0, Upper Inf), data dat)关键参数是Lower和Upper它们定义了审查边界注意是首字母大写。这一点和AER::tobit()完全不同那个函数用的是小写的left和right。我见过太多人把两套参数混着用结果模型静默跑完系数却完全不是一回事。设定规则如下左审查Lower 0表示所有 y ≤ 0 的观测被审查为0右审查Upper 10000表示所有 y ≥ 10000 的观测被审查为10000双审查Lower 0, Upper 10000同时设置单侧无审查另一侧填Inf或-Inf。注意VGAM 不需要额外提供一个“是否被审查”的哑变量。你只要把实际观测到的 y 传给公式审查点上的重复取值就是信息的一部分它会自动进入似然函数。2.3 σ的第二条线性预测器与zero参数普通回归只需要估计一组系数 β。tobit 还需要估计 σ因为审查概率的表达式里含有 σ。VGAM 把 σ 作为第二条线性预测器linear predictor来处理默认使用 log 链接函数所以估计的是 log(σ) 而不是 σ。这就是为什么summary(fit)或者coef(fit)里会出现一个(Intercept):2这样的项它的含义是截距对 log(σ) 的贡献拟合后需要做指数变换才能得到 σ。默认的tobit()族只让 σ 是常数也就是同方差假设。如果你怀疑审查机制和异方差有关可以尝试调整 family 里的zero参数。zero在 VGAM 里用来指定哪条线性预测器只保留截距项。把对应 σ 的预测器从“只含截距”改为“允许协变量影响”就变成一个异方差tobit模型。这种做法要注意模型会复杂很多对数据量的要求也更高不建议一上来就放开。2.4 应变量单位和收敛容差tobit 的似然函数是 β 和 σ 联合估计的σ 通过 log 链接进入优化过程。如果应变量量级特别大比如原始单位是“分”数值动辄几百万log(σ) 的初始值很可能离最优解很远导致迭代不稳定。实操中我倾向先把应变量统一到一个好读的单位。比如金额用“万元”而不是“元”理赔额用“千元”。这样估计出来的 β 虽然数值会变但模型解释、边际效应都不会受到影响。单位没有调整清楚之前先不要急着怪 VGAM 不收敛。3. 从模拟数据开始跑通一次vglm拟合3.1 生成一个有真实答案的审查数据为了看清楚 VGAM 的估计效果我们造一份已知真实系数的数据。设潜变量模型为y* 0.8 1.2·x1 - 0.6·x2 εε ~ N(0, 2²)然后左审查在0处观测到的 y max(0, y*)。set.seed(20240915) n - 5000 x1 - rnorm(n) x2 - rnorm(n, 2, 1) ystar - 0.8 1.2 * x1 - 0.6 * x2 rnorm(n, 0, 2) y - ifelse(ystar 0, ystar, 0) dat - data.frame(x1, x2, y) cat(审查比例, mean(y 0), \n)现实中审查比例没有硬性上限但审查比例太高时信息量会明显下降。这里的模拟大概会有两到三成的0值比较接近常见调查数据的形态。3.2 用vglm拟合tobit并解读summarylibrary(VGAM) fit - vglm(y ~ x1 x2, tobit(Lower 0, Upper Inf), data dat) summary(fit)输出里最核心的是系数表。因为这个 family 有两条线性预测器你会看到两类系数(Intercept):1、x1、x2潜变量方程 y* Xβ 里的截距和系数(Intercept):2log(σ) 的截距。我们模拟时设置的真实值是 β0 0.8β1 1.2β2 -0.6σ 2。拟合出来之后x1的系数应该接近1.2x2接近-0.6(Intercept):1接近0.8。而 σ 要通过exp(coef(fit)[(Intercept):2])转回来结果应该接近2。手动提取时建议这样写b - coef(fit) sigma_hat - exp(b[[(Intercept):2]]) print(sigma_hat)如果直接在摘要里找 σ很多第一次用的人会误把(Intercept):2当成另一个截距这不对。3.3 和AER::tobit的结果做对照VGAM 的结果不一定能让你完全安心常见做法是再用AER::tobit()跑一遍两相对照。library(AER) fit_aer - tobit(y ~ x1 x2, left 0, right Inf, data dat) summary(fit_aer)结果出来后系数部分和vglm()基本一致差异一般出现在小数后第三四位属于优化算法和收敛容差的正常波动。AER::tobit()输出的Log(scale)就是 log(σ)和 VGAM 的(Intercept):2含义相同。我有一个习惯正式出结果前用两个包互相验证。如果两边系数差别很大第一件事不是查统计学理论而是检查审查参数是不是写反了尤其是Lower/Upper和left/right混用。3.4 三种边际效应潜变量、观测值、审查概率tobit 的系数不能像 OLS 那样直接解读成“x 增加一单位y 平均变化多少”。因为真正的效应分三种第一对潜变量 y* 的条件均值 E(y*|x) 的效应等于 βj。第二对观测值 y 的条件均值 E(y|x) 的效应。在左审查点为0的情况下它的公式可以化简为∂E(y|x)/∂xj βj · Φ(xβ/σ)也就是系数乘上“潜变量超过审查点的概率”。这比 βj 小因为审查机制把一部分效应吸收掉了。第三对“y 0 的概率”的效应∂P(y0|x)/∂xj (βj/σ) · φ(xβ/σ)很多报告只写第一类这是不完整的。如果业务关心的是实际观测支出应该报告第二类。如果关心的是“是否购买”应该报告第三类。在 R 里可以自己写一个小函数用均值处的 x 来计算b - coef(fit) beta - c(b[[x1]], b[[x2]]) sigma - exp(b[[(Intercept):2]]) xb - b[[(Intercept):1]] b[[x1]] * dat$x1 b[[x2]] * dat$x2 Phi - pnorm(xb / sigma) phi - dnorm(xb / sigma) me_obs - beta * mean(Phi) me_prob - beta / sigma * mean(phi)这样输出的me_obs是观测支出均值的平均边际效应me_prob是“支出为正概率”的平均边际效应。3.5 计算预测期望值vglm()的predict()默认返回的是线性预测子不是观测支出 y 的期望。很多人在这一步踩坑拿到的预测值全是负数怀疑模型坏了。线性预测子的第一列是 μ xβ第二列是 log(σ)。要得到观测支出 y 的期望需要手动套公式E(y|x) Φ(xβ/σ) · xβ σ · φ(xβ/σ)代码可以这么写lp - predict(fit, newdata dat) mu - lp[, 1] sigma_pred - exp(lp[, 2]) Phi - pnorm(mu / sigma_pred) phi - dnorm(mu / sigma_pred) yhat_obs - Phi * mu sigma_pred * phi这样得到的yhat_obs才是“观测到的 y”的模型预测值。和实际数据对比时要用这个值而不是predict()的第一列。4. 怎么看拟合结果诊断和模型比较4.1 对数似然、AIC与似然比检验tobit 是用极大似然估计的所以模型比较的基础也是似然。VGAM 里可以直接用AIC(fit)、BIC(fit)来比较非嵌套模型。嵌套模型比较最常用的是似然比检验。比如想判断 x2 是否需要进入模型fit_small - vglm(y ~ x1, tobit(Lower 0, Upper Inf), data dat) lmtest::lrtest(fit_small, fit)如果检验统计量显著说明加入 x2 显著提升了模型拟合。注意tobit 模型的 AIC/BIC 计算使用了完整似然包括审查点的累积概率部分所以不要只拿正数部分的拟合优度去比。4.2 残差不能像普通回归那样看普通线性回归里残差服从正态分布、随机散布是模型正确的标志。tobit 模型里所有 y 0 的观测残差都是负数集中堆积这种堆叠是审查机制内生的不代表模型不好。我建议的诊断思路有两种。第一把 y 0 的观测单独拿出来计算标准化残差并画出 QQ 图检查正态性假设是否合理。审查点以下的观测无法计算连续残差但它们的信息已经通过累积概率进入了似然。第二比较模型预测的审查比例与实际审查比例。用拟合参数计算 P(y 0|x) 的平均值看看是否接近样本里的实际非零比例。如果差很多说明模型对审查概率的刻画有问题。4.3 用模拟数据检验恢复效果模拟数据的最大优势是我们可以拿估计值和真实值对表。以刚才的数据为例如果反复重复几百次模拟估计值的均值应该接近真实参数估计值的标准差应该接近 summary 里的标准误。这比任何残差图都更能说明模型的正确性。实际工作中如果数据不是模拟出来的没有真实参数可对照可以在分析前人为抽取一个子集把某些观测的 y 改成审查点重新拟合看看参数变化是否在合理范围内。这类敏感性分析对受限因变量模型特别有用因为审查比例直接影响估计精度。5. 真实数据分析中的几道坎查错清单5.1 审查方向写反了症状与自检用 VGAM 跑 tobit最隐蔽、最容易出问题的就是审查方向。Lower 0表示左审查意味着把低于0的真实值记录为0。如果你的数据是“支出金额”左审查在0处这是对的。如果数据是“最高收入上限”比如超出50000元统一记录为50000那应该设置Upper 50000同时Lower -Inf或设为Lower -Inf。写反之后模型仍会正常迭代但系数会发生系统性偏移符号甚至可能反向。自查方法很简单计算 y 等于审查边界的比例。如果左审查设了0但数据根本没有0只在50000处堆积那一定是审查方向不对。另外把 VGAM 的结果和AER::tobit()对照是更快的检查方式。5.2 不收敛缩放、迭代次数、初始值vglm()偶尔会报“Iterations terminated because half-step sizes were very small”之类的警告。这通常不是模型理论问题而是数值优化问题。我处理过的情况最常见的原因是应变量量级太大。单位从“元”换成“万元”问题立刻消失。第二种原因是协变量量级差异过大比如一个变量在0到1之间另一个变量在几万量级模型矩阵条件数很差。先标准化连续变量再重新拟合通常能解决。如果仍然不收敛可以调大迭代次数fit - vglm(y ~ x1 x2, tobit(Lower 0, Upper Inf), data dat, maxit 100, trace TRUE)trace TRUE可以在迭代过程中输出对数似然方便观察是否在缓慢爬升。合适的初始值也很重要但 VGAM 自带的初始化对大多数问题已经足够不建议手动乱给。5.3 审查点数值不能随便改有些数据整理时会把审查点替换成其他值。比如左审查点明明是0有人为了“方便计算”把0替换成0.01或者把右审查点替换成上限的1.1倍。这样做等于给模型注入了虚假信息估计结果会产生系统性偏移。正确做法是原始值是0就保留0原始值是上限就保留上限。tobit 模型不需要你“修正”审查数据它要的就是审查点上的堆积信息。5.4 异方差tobit与扩展思路经典 tobit 假设 ε 的方差是常数。如果实际数据的离散程度随 x 增大而增大σ 就不是常数此时简单 tobit 的 β 可能仍然是近似一致的但标准误会偏推断结论不可靠。VGAM 的优势在于可以扩展这种设定通过调整zero参数让 log(σ) 那条线性预测器也包含协变量。不过接下来要注意模型不再只有一个 σ而是一组随协变量变化的 σ(x)。这对数据量的要求高很多解释起来也更复杂。我的建议是先用简单同方差模型把基准结果做出来一旦发现审查比例在各 x 分组下差异很大再去尝试异方差设定。不要一上来就把模型复杂度拉满否则容易被优化问题和解释困难拖住。还有一个容易混淆的边界问题tobit 假设“是否被审查”和“观测值大小”由同一个潜变量机制决定。如果现实里的机制是两阶段比如先决定“是否就业”再决定“就业后的工资水平”这两个决策的影响因素可能是不同的。这时候应该考虑样本选择模型而不是单纯套 tobit。VGAM 虽然功能强大但也不是万能的模型选择永远要服从业务机制。最后再分享一条实际操作经验每次拟合完 tobit我都会同时保留三样东西——拟合参数的原始输出、预测的潜变量均值、预测的观测值期望。业务报告里通常需要观测值期望而方法学验证需要潜变量均值只保留其中一个后面要做边际效应或绘图时又会重新跑一遍模型白白浪费时间。