我的问题是关于不必要的预测器,即不提供任何新的线性信息的变量,或者是其他预测器的线性组合的变量。如您所见,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
。
有人能解释一下原因吗?
发布于 2017-06-23 04:27:34
会有
NA
吗?
是。添加这些列不会扩大列空间。得到的矩阵是秩亏的.
多少
NA
?
这取决于数字的等级。
number of NA = number of coefficients - rank of model matrix
在您的示例中,在引入ec
之后,将有一个NA
。更改模型公式中协变量的规格顺序实质上是对模型矩阵进行列改组。这不会改变矩阵的排序,因此无论您的规范顺序如何,您总是只得到一个NA
。
好的,但是哪一个是
NA
?
lm
使用限制的列旋转进行LINPACK QR分解。协变量的顺序影响到哪一个是NA
。一般来说,“先到先得”的原则是成立的,NA
的地位是相当可预测的。以你的例子为例。在第一个规范中,这些协线性项以Examination
,Catholic
,ec
顺序出现,因此第三个ec
具有NA
系数。在第二个规范中,这些术语以ec
、Examination
、Catholic
顺序显示,第三个Catholic
具有NA
系数。请注意,系数估计并不是不变量的规格顺序,虽然拟合的值是不变的。
如果采用完全列旋转的LAPACK分解,则系数估计将不受规范阶的影响。然而,NA
的位置并不像在LINPACK情况下那样可预测,它完全是由数字决定的。
数值算例
在mgcv
包中实现了基于LAPACK的QR分解。当使用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
和bam
估计下的位置。
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
的看法相同。
发布于 2017-06-23 04:36:29
ec、检测和是3个参数,需要至少2个变量来确定第三个。重要的是,每三个中有2个是必需的。现在,当你把这个传递给lm时,3个相关变量中的前两个将得到系数,第三个变量将以NA结束。变量的顺序很重要。我希望这解释了考试和天主教两者都不是NA的原因。只有电子商务,你不能同时决定考试和天主教。
https://stackoverflow.com/questions/44721341
复制