由于具有不必要变量的线性回归，R 是否总是 return NA 作为系数？

Question

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

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

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

ec <- swiss$Examination + swiss$Catholic

当我们运行一个包含不必要变量的线性回归时，R 会删除其他项的线性组合项和 returns 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.

谁能解释一下这是什么原因？

Answer 1

There will be NA?

是的。添加这些列不会扩大列 space。生成的矩阵秩亏。

How many NA?

这取决于数字排名。

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

在你的例子中，引入ec后，就会有一个NA。更改模型公式中协变量的指定顺序实质上是对模型矩阵进行列改组。这不会改变矩阵秩，因此无论您的指定顺序如何，您总是只会得到一个 NA。

OK, but which one is NA?

lm 使用 restricted 列旋转进行 LINPACK QR 分解。协变量的顺序影响哪一个是 NA。一般来说，“先到先得”的原则是成立的，NA的位置是可以预见的。以你的例子为例。在第一个规范中，这些共线性项以 Examination、Catholic、ec 的顺序出现，因此第三个 ec 具有 NA 系数。在您的第二个规范中，这些项以 ec、Examination、Catholic 顺序显示，第三个 Catholic 具有 NA 系数。请注意，尽管拟合值是不变的，但系数估计并非对规范顺序不变。

如果采用 LAPACK QR 因式分解 complete 列旋转，系数估计将不随规范顺序变化。然而，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 的列，看看 NA 是否在 lm 估计下改变了它的位置，或者 0 是否在 gam 和 [=40] 下改变了它的位置=]估计。

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)

我们看到 NA 通过 X 的列改组改变了它的位置。估计系数也不同。

现在我们检查 gam

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

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

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

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

结果与我们看到的 gam.

相同

Answer 2

ec , examination 和 catholic 是您需要的 3 个参数 至少有2个变量来确定第三个。重要的是总是需要三分之二。现在，当您将其传递给 lm 时，3 个相关变量中的前两个将获得系数，第三个将以 NA 结束。变量的顺序很重要。我希望这可以解释为什么 examination 和 catholic 都不适用。仅凭 ec，您无法同时确定考试和天主教

由于具有不必要变量的线性回归，R 是否总是 return NA 作为系数？

Does R always return NA as a coefficient as a result of linear regression with unnecessary variables?

regression

r

linear-regression

lm

coefficients

数值示例