R语言中poly()拟合多项式回归后调用vif()报错的问题咨询

问题场景 #

当使用poly()函数在R中拟合多元二次多项式回归，随后调用vif()计算方差膨胀因子时，会触发如下错误：

y = c(0.22200,0.39500,0.42200,0.43700,0.42800,0.46700,0.44400,0.37800,0.49400,
      0.45600,0.45200,0.11200,0.43200,0.10100,0.23200,0.30600,0.09230,0.11600,
      0.07640,0.43900,0.09440,0.11700,0.07260,0.04120,0.25100,0.00002)
x1 = c(7.3,8.7,8.8,8.1,9.0,8.7,9.3,7.6,10.0,8.4,9.3,7.7,9.8,7.3,8.5,9.5,7.4,7.8,
       7.7,10.3,7.8,7.1,7.7,7.4,7.3,7.6)
x2 = c(0.0,0.0,0.7,4.0,0.5,1.5,2.1,5.1,0.0,3.7,3.6,2.8,4.2,2.5,2.0,2.5,2.8,2.8,
       3.0,1.7,3.3,3.9,4.3,6.0,2.0,7.8)
x3 = c(0.0,0.3,1.0,0.2,1.0,2.8,1.0,3.4,0.3,4.1,2.0,7.1,2.0,6.8,6.6,5.0,7.8,7.7,
       8.0,4.2,8.5,6.6,9.5,10.9,5.2,20.7)
m = lm(y~poly(x1, x2, x3, degree=2, raw=TRUE))
summary(m)

调用vif()时的错误：

> vif(m)
Error in vif.default(m) : model contains fewer than 2 terms

但实际模型包含9个自变量项+1个截距，m$rank返回结果为10，证实模型确实有多个项。

问题原因 #

你的怀疑部分正确，但本质不是poly()与vif()不兼容，而是**poly()在lm公式中被当作单一的矩阵型预测变量**，而非多个独立的自变量。vif()函数会将整个poly()输出的矩阵识别为一个变量，因此判断模型自变量数不足2，抛出错误。

解决办法 #

方法1：手动展开多项式项

直接在公式中写出所有二次项（主效应、平方项、交互项），让每个项成为独立变量：

m2 = lm(y ~ x1 + x2 + x3 + I(x1^2) + I(x2^2) + I(x3^2) + I(x1*x2) + I(x1*x3) + I(x2*x3))
vif(m2)

这样vif()能识别每个独立的自变量项，正常计算VIF值。

方法2：拆分poly()输出为单独变量

先将poly()生成的矩阵转换为数据框，合并到原数据中，再拟合模型：

# 生成多项式矩阵并命名
poly_terms = as.data.frame(poly(x1, x2, x3, degree=2, raw=TRUE))
colnames(poly_terms) = paste0("term_", 1:ncol(poly_terms))

# 合并数据
data_all = cbind(data.frame(y, x1, x2, x3), poly_terms)

# 拟合模型（去掉原x1/x2/x3避免重复）
m3 = lm(y ~ ., data = data_all[, !colnames(data_all) %in% c("x1", "x2", "x3")])
vif(m3)

方法3：改进手动计算方式

你给出的手动计算代码存在偏差，因为VIF的定义是1/(1-Rj²)（Rj²是第j个自变量对其余所有自变量回归的决定系数），而非直接取(X'X)逆矩阵对角线的倒数。正确的手动计算方式如下：

# 获取不含截距的自变量矩阵
X = model.matrix(m)[, -1]
vifs = numeric(ncol(X))

for (i in 1:ncol(X)) {
  # 将第i个变量作为因变量，对其余变量回归
  r_sq = summary(lm(X[, i] ~ X[, -i]))$r.squared
  vifs[i] = 1 / (1 - r_sq)
}
names(vifs) = colnames(X)
vifs

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