对于具有线性依赖关系的模型,Reqsubsets 结果与 coef() 不同
Reqsubsets results differ with coef() for model with linear dependencies
在对具有线性相关性的数据使用包 leaps
中的 Regsubsets
时,我发现 coef()
和 summary()$which
给出的结果不同。似乎,当找到线性相关性时,重新排序会改变系数的位置和 coef()
returns 错误的值。
我使用 mtcars
只是为了 "simulate" 我在处理其他数据时遇到的问题。在第一个例子中没有 lin 的问题。 BIC 的依赖项和最佳给定模型是 mpg~wt+cyl
并且 coef()
、summary()$which
都给出相同的结果。在第二个示例中,我添加了虚拟变量,因此存在完全多重共线性的可能性,但按此顺序排列的变量(最后一列中的虚拟变量)不会导致问题。在最后一个示例中,在更改数据集中变量的顺序后,问题终于出现了,coef()
、summary()$which
给出了不同的模型。这种方法有什么不正确的地方吗?还有其他方法可以从 regsubsets
获取系数吗?
require("leaps") #install.packages("leaps")
###Example1
dta <- mtcars[,c("mpg","cyl","am","wt","hp") ]
bestSubset.cars <- regsubsets(mpg~., data=dta)
(best.sum <- summary(bestSubset.cars))
#
w <- which.min(best.sum$bic)
best.sum$which[w,]
#
best.sum$outmat
coef(bestSubset.cars, w)
#
###Example2
dta2 <- cbind(dta, manual=as.numeric(!dta$am))
bestSubset.cars2 <- regsubsets(mpg~., data=dta)
(best.sum2 <- summary(bestSubset.cars2))
#
w <- which.min(best.sum2$bic)
best.sum2$which[w,]
#
coef(bestSubset.cars2, w)
#
###Example3
bestSubset.cars3 <- regsubsets(mpg~., data=dta2[,c("mpg","manual","am","cyl","wt","hp")])
(best.sum3 <- summary(bestSubset.cars3))
#
w <- which.min(best.sum3$bic)
best.sum3$which[w,]
#
coef(bestSubset.cars3, w)
#
best.sum2$which
coef(bestSubset.cars2,1:4)
best.sum3$which
coef(bestSubset.cars3,1:4)
summary.regsubsets 和 regsubsets 的 vars 顺序不同。 regsubsets 的通用函数 coef() 在一个函数中调用这两个函数,如果您尝试 force.in 或使用固定顺序的公式,结果会很混乱。更改 coef() 函数中的某些行可能会有所帮助。试试下面的代码,看看它是否有效!
coef.regsubsets <- function (object, id, vcov = FALSE, ...)
{
s <- summary(object)
invars <- s$which[id, , drop = FALSE]
betas <- vector("list", length(id))
for (i in 1:length(id)) {
# added
var.name <- names(which(invars[i, ]))
thismodel <- which(object$xnames %in% var.name)
names(thismodel) <- var.name
# deleted
#thismodel <- which(invars[i, ])
qr <- .Fortran("REORDR", np = as.integer(object$np),
nrbar = as.integer(object$nrbar), vorder = as.integer(object$vorder),
d = as.double(object$d), rbar = as.double(object$rbar),
thetab = as.double(object$thetab), rss = as.double(object$rss),
tol = as.double(object$tol), list = as.integer(thismodel),
n = as.integer(length(thismodel)), pos1 = 1L, ier = integer(1))
beta <- .Fortran("REGCF", np = as.integer(qr$np), nrbar = as.integer(qr$nrbar),
d = as.double(qr$d), rbar = as.double(qr$rbar), thetab = as.double(qr$thetab),
tol = as.double(qr$tol), beta = numeric(length(thismodel)),
nreq = as.integer(length(thismodel)), ier = numeric(1))$beta
names(beta) <- object$xnames[qr$vorder[1:qr$n]]
reorder <- order(qr$vorder[1:qr$n])
beta <- beta[reorder]
if (vcov) {
p <- length(thismodel)
R <- diag(qr$np)
R[row(R) > col(R)] <- qr$rbar
R <- t(R)
R <- sqrt(qr$d) * R
R <- R[1:p, 1:p, drop = FALSE]
R <- chol2inv(R)
dimnames(R) <- list(object$xnames[qr$vorder[1:p]],
object$xnames[qr$vorder[1:p]])
V <- R * s$rss[id[i]]/(object$nn - p)
V <- V[reorder, reorder]
attr(beta, "vcov") <- V
}
betas[[i]] <- beta
}
if (length(id) == 1)
beta
else betas
}
另一个对我有用的解决方案是在 运行 regsubsets 之前随机化数据集中列(自变量)的顺序。这个想法是,在重新排序后,希望高度相关的列彼此相距很远,并且不会触发 regsubsets 算法中的重新排序行为。
在对具有线性相关性的数据使用包 leaps
中的 Regsubsets
时,我发现 coef()
和 summary()$which
给出的结果不同。似乎,当找到线性相关性时,重新排序会改变系数的位置和 coef()
returns 错误的值。
我使用 mtcars
只是为了 "simulate" 我在处理其他数据时遇到的问题。在第一个例子中没有 lin 的问题。 BIC 的依赖项和最佳给定模型是 mpg~wt+cyl
并且 coef()
、summary()$which
都给出相同的结果。在第二个示例中,我添加了虚拟变量,因此存在完全多重共线性的可能性,但按此顺序排列的变量(最后一列中的虚拟变量)不会导致问题。在最后一个示例中,在更改数据集中变量的顺序后,问题终于出现了,coef()
、summary()$which
给出了不同的模型。这种方法有什么不正确的地方吗?还有其他方法可以从 regsubsets
获取系数吗?
require("leaps") #install.packages("leaps")
###Example1
dta <- mtcars[,c("mpg","cyl","am","wt","hp") ]
bestSubset.cars <- regsubsets(mpg~., data=dta)
(best.sum <- summary(bestSubset.cars))
#
w <- which.min(best.sum$bic)
best.sum$which[w,]
#
best.sum$outmat
coef(bestSubset.cars, w)
#
###Example2
dta2 <- cbind(dta, manual=as.numeric(!dta$am))
bestSubset.cars2 <- regsubsets(mpg~., data=dta)
(best.sum2 <- summary(bestSubset.cars2))
#
w <- which.min(best.sum2$bic)
best.sum2$which[w,]
#
coef(bestSubset.cars2, w)
#
###Example3
bestSubset.cars3 <- regsubsets(mpg~., data=dta2[,c("mpg","manual","am","cyl","wt","hp")])
(best.sum3 <- summary(bestSubset.cars3))
#
w <- which.min(best.sum3$bic)
best.sum3$which[w,]
#
coef(bestSubset.cars3, w)
#
best.sum2$which
coef(bestSubset.cars2,1:4)
best.sum3$which
coef(bestSubset.cars3,1:4)
summary.regsubsets 和 regsubsets 的 vars 顺序不同。 regsubsets 的通用函数 coef() 在一个函数中调用这两个函数,如果您尝试 force.in 或使用固定顺序的公式,结果会很混乱。更改 coef() 函数中的某些行可能会有所帮助。试试下面的代码,看看它是否有效!
coef.regsubsets <- function (object, id, vcov = FALSE, ...)
{
s <- summary(object)
invars <- s$which[id, , drop = FALSE]
betas <- vector("list", length(id))
for (i in 1:length(id)) {
# added
var.name <- names(which(invars[i, ]))
thismodel <- which(object$xnames %in% var.name)
names(thismodel) <- var.name
# deleted
#thismodel <- which(invars[i, ])
qr <- .Fortran("REORDR", np = as.integer(object$np),
nrbar = as.integer(object$nrbar), vorder = as.integer(object$vorder),
d = as.double(object$d), rbar = as.double(object$rbar),
thetab = as.double(object$thetab), rss = as.double(object$rss),
tol = as.double(object$tol), list = as.integer(thismodel),
n = as.integer(length(thismodel)), pos1 = 1L, ier = integer(1))
beta <- .Fortran("REGCF", np = as.integer(qr$np), nrbar = as.integer(qr$nrbar),
d = as.double(qr$d), rbar = as.double(qr$rbar), thetab = as.double(qr$thetab),
tol = as.double(qr$tol), beta = numeric(length(thismodel)),
nreq = as.integer(length(thismodel)), ier = numeric(1))$beta
names(beta) <- object$xnames[qr$vorder[1:qr$n]]
reorder <- order(qr$vorder[1:qr$n])
beta <- beta[reorder]
if (vcov) {
p <- length(thismodel)
R <- diag(qr$np)
R[row(R) > col(R)] <- qr$rbar
R <- t(R)
R <- sqrt(qr$d) * R
R <- R[1:p, 1:p, drop = FALSE]
R <- chol2inv(R)
dimnames(R) <- list(object$xnames[qr$vorder[1:p]],
object$xnames[qr$vorder[1:p]])
V <- R * s$rss[id[i]]/(object$nn - p)
V <- V[reorder, reorder]
attr(beta, "vcov") <- V
}
betas[[i]] <- beta
}
if (length(id) == 1)
beta
else betas
}
另一个对我有用的解决方案是在 运行 regsubsets 之前随机化数据集中列(自变量)的顺序。这个想法是,在重新排序后,希望高度相关的列彼此相距很远,并且不会触发 regsubsets 算法中的重新排序行为。