企业官网建设流程全解析

之前有个项目要处理一份个人消费问卷，样本量三千多，因变量是“过去三个月网购花费”。数据里四成左右的人填了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并解读summary

library(VGAM) fit <- vglm(y ~ x1 + x2, tobit(Lower = 0, Upper = Inf), data = dat) summary(fit)

输出里最核心的是系数表。因为这个 family 有两条线性预测器，你会看到两类系数：

• (Intercept):1、x1、x2：潜变量方程 y* = Xβ 里的截距和系数；
• (Intercept):2：log(σ) 的截距。

我们模拟时设置的真实值是 β0 = 0.8，β1 = 1.2，β2 = -0.6，σ = 2。拟合出来之后，x1的系数应该接近1.2，x2接近-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(y>0|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，我都会同时保留三样东西——拟合参数的原始输出、预测的潜变量均值、预测的观测值期望。业务报告里通常需要观测值期望，而方法学验证需要潜变量均值，只保留其中一个，后面要做边际效应或绘图时又会重新跑一遍模型，白白浪费时间。