如何在R中运行Box-Tidwell检验验证逻辑回归自变量与对数几率的线性关系

Box-Tidwell检验的核心作用是验证连续自变量与模型log odds的线性关系假设，检验仅针对连续型变量开展，分类变量BURN不需要纳入检验流程。

操作前准备 #

• 完成数据集基础清洗，剔除缺失值、明显异常值
• 确保分类变量BURN在R中被识别为因子类型，避免被误判为连续变量参与计算
• 提前整理好5个连续自变量：林分平均年龄AGE、林冠覆盖率CANCOV、地被覆盖率GROCOV、林分面积AREA、林分断面积BASAR

分步操作流程 #

1. 拟合基线逻辑回归模型 #

先不加入任何检验交互项，用所有预测变量搭建基础模型：

# 将your_dataset替换为你的实际数据集名称，将“啄木鸟存在缺失”替换为你实际记录存在/缺失的二分类响应变量名
base_model <- glm(
  formula = 啄木鸟存在缺失 ~ AGE + CANCOV + GROCOV + AREA + BASAR + factor(BURN),
  family = binomial(link = "logit"),
  data = your_dataset
)

注：公式中对BURN使用factor()包裹，是为了明确告知R这是三分类变量，不会把0/1/2的取值识别为连续梯度参与模型计算

2. 构造检验所需的交互项 #

Box-Tidwell检验的核心逻辑是给每个连续自变量，构造「变量值 * 该变量值的自然对数」形式的新交互项。如果这些交互项的系数不显著，就说明原变量和log odds的线性关系成立。

# 为每个连续变量生成对应交互项，加1e-6极小偏移避免log(0)报错
bt_test_data <- your_dataset %>%
  mutate(
    AGE_inter = AGE * log(AGE + 1e-6),
    CANCOV_inter = CANCOV * log(CANCOV + 1e-6),
    GROCOV_inter = GROCOV * log(GROCOV + 1e-6),
    AREA_inter = AREA * log(AREA + 1e-6),
    BASAR_inter = BASAR * log(BASAR + 1e-6)
  )

注：加1e-6的偏移量是为了规避变量取值为0时，自然对数计算出现无穷大的报错，该偏移量极小，不会对检验结果的判断造成干扰

3. 拟合检验模型并提取结果 #

把构造好的5个交互项全部加入基线模型，输出结果后重点查看交互项的显著性：

# 拟合带所有交互项的扩展检验模型
bt_test_model <- glm(
  formula = 啄木鸟存在缺失 ~ AGE + CANCOV + GROCOV + AREA + BASAR + factor(BURN) +
    AGE_inter + CANCOV_inter + GROCOV_inter + AREA_inter + BASAR_inter,
  family = binomial(link = "logit"),
  data = bt_test_data
)

# 输出完整模型结果
summary(bt_test_model)

结果判断规则 #

• 逐一查看模型输出中后缀为_inter的交互项对应的p值
• 若某交互项p值 > 0.05，说明对应的连续自变量满足与log odds的线性关系假设，不需要对该变量做额外转换
• 若某交互项p值 < 0.05，说明对应的连续自变量与log odds存在非线性关系，需要对该变量做对数转换、分箱处理、加入二次多项式项等操作修正模型设定
• 可配合似然比检验对比基线模型和带交互项模型的解释度差异，辅助判断整体线性假设是否成立

实操注意事项 #

• 不要将分类变量BURN纳入Box-Tidwell检验流程，分类变量本身不满足连续尺度的线性假设检验前提
• 当变量值跨度较大（比如林分面积AREA从0.1公顷到上千公顷）时，检验出非线性关系的概率会升高，不要机械套用线性假设
• 当交互项p值处于0.05的临界值附近时，可搭配loess局部拟合曲线、分变量的残差图辅助判断，不要完全依赖p值下结论

内容的提问来源于stack exchange，提问作者Casey