在进行线性回归时，如果存在不必要的变量，是否始终会得到R作为系数的NA值？

Question

我的问题是关于不必要的预测变量，也就是那些不提供任何新的线性信息或者是其他预测变量的线性组合。正如你所见，swiss 数据集有六个变量。

library(swiss)
names(swiss)
# "Fertility"        "Agriculture"      "Examination"      "Education"        
# "Catholic"      "Infant.Mortality"

现在我介绍一个新变量ec。它是Examination和Education的线性组合。

ec <- swiss$Examination + swiss$Catholic

当我们运行不必要变量的线性回归时，R会删除线性组合其他项的项，并将其系数返回为NA。下面的命令完美地说明了这一点。

lm(Fertility ~ . + ec, swiss)

Coefficients:
 (Intercept)       Agriculture       Examination         Education            
     66.9152           -0.1721           -0.2580           -0.8709 

Catholic  Infant.Mortality    ec

  0.1041            1.0770    NA

然而，当我们按照下面所示的方法先回归ec，然后再回归所有的自变量时，

lm(Fertility ~ ec + ., swiss)

 Coefficients:
 (Intercept)                ec       Agriculture       Examination           
     66.9152            0.1041           -0.1721           -0.3621           
  Education          Catholic     Infant.Mortality  
    -0.8709                NA            1.0770  

我希望两个系数Catholic和Examination都应该是NA。变量ec是它们的线性组合，但最终Examination的系数不是NA，而Catholic的系数是NA。

有人可以解释一下原因吗？

- BRCN

---

如果您有三个变量A，B，A+B（其中A和B本身不共线），则可以估计任意两个变量的参数。 - Ben Bolker

这是线性代数的基本原则。您有一个排名为2的矩阵x1、x2，现在您添加一个变量x3，它是x1和x2的线性组合。该矩阵仍然具有排名2。 - lmo

Answer 1

会有 NA 吗？

是的。添加这些列不会增加列空间，导致矩阵的秩降低。

有多少个 NA？

这取决于数值秩。

number of NA = number of coefficients - rank of model matrix

在你的示例中，引入ec后，将会有一个NA。在模型公式中更改协变量的规范顺序实质上是对模型矩阵进行列重排。这不会改变矩阵秩，因此无论您的规范顺序如何，始终只会得到一个NA。 lm使用受限列置换的LINPACK QR分解。协变量的顺序影响哪个是NA。一般地，“先来先服务”的原则适用，并且NA的位置非常可预测。以你的示例为例，在第一种规范下，这些共线项按照Examination，Catholic，ec的顺序出现，因此第三个ec具有NA系数。在第二个规范中，这些项按照ec，Examination，Catholic的顺序出现，因此第三个Catholic具有NA系数。请注意，系数估计与规范顺序不变，尽管拟合值是不变的。

如果采用完全列置换的LAPACK QR分解，则系数估计将相对于规范顺序不变。然而，NA的位置不如在LINPACK情况下可预测，并且完全由数值决定。

---

数值示例

基于LAPACK的QR分解实现在mgcv包中。当使用REML估计时检测数值秩，并将不可识别的系数报告为0（而不是NA）。因此，我们可以比较线性模型估计中的lm和gam/bam。让我们首先构造一个玩具数据集。

set.seed(0)

# an initial full rank matrix
X <- matrix(runif(500 * 10), 500)
# make the last column as a random linear combination of previous 9 columns
X[, 10] <- X[, -10] %*% runif(9)

# a random response
Y <- rnorm(500)

现在我们对 X 的列进行洗牌，以查看在 lm 估计下 NA 是否改变其位置，或者在 gam 和 bam 估计下 0 是否改变其位置。

test <- function (fun = lm, seed = 0, ...) {
  shuffleFit <- function (fun) {
    shuffle <- sample.int(ncol(X))
    Xs <- X[, shuffle]
    b <- unname(coef(fun(Y ~ Xs, ...)))
    back <- order(shuffle)
    c(b[1], b[-1][back])
    }
  set.seed(seed)
  oo <- t(replicate(10, shuffleFit(fun)))
  colnames(oo) <- c("intercept", paste0("X", 1:ncol(X)))
  oo
  }

首先，我们使用 lm 进行检查：

test(fun = lm)

我们发现，在对 X 进行列重排时，NA 的位置也会改变。估计系数也会有所变化。

---

现在我们使用 gam 进行检查。

library(mgcv)
test(fun = gam, method = "REML")

我们看到估计值对于 X 的列置换是不变的，X5 的系数始终为0。

---

最后我们检查 bam（bam 对于像这样的小数据集来说速度较慢。它设计用于大型或超大型数据集。所以下面的速度明显较慢）。

test(fun = bam, gc.level = -1)

这个结果与我们在gam中看到的一样。

- Zheyuan Li

Answer 2

ec、examination和catholic是三个参数，你需要至少2个变量来确定第三个变量。重要的是，这三个参数中必须至少有两个变量。当你将这些参数传递给lm时，前两个相关变量将获得系数，第三个变量将以NA结尾。变量的顺序非常重要。希望这能解释为什么考试和天主教都不是NA。只有ec一个变量是无法同时确定考试和天主教的。

- niths4u