将指示符矩阵对应于 R 中的列名

Question

我要解决的问题是“如何创建一系列自动化代码，从数据集中提取所需的列标题名称以适应一般线性化模型 (glm)？”我有一个包含 8 个变量的数据集；但是，我只想使用其中的 3 个来交叉验证并找到“最佳”模型。这是我想出的：

library(boot)

salary <- read.csv("salary_data.csv")
vars <- colnames(salary[c(2,3,7)])
nvars <- length(vars)
list.to.expand = vector(mode = "list", length = nvars)
for (i in 1:nvars){
  list.to.expand[[i]]=c(0,1)
}
model.spec.matrix <- expand.grid(list.to.expand)

vars
model.spec.matrix

names(model.spec.matrix) <- vars
vars.to.use <- model.spec.matrix[2,]
vars.to.use <- as.numeric(vars.to.use)
cn <- c()
for (i in 1:nrow(model.spec.matrix)){
  if(i==1){cn <- colnames(model.spec.matrix[sapply(model.spec.matrix[i,], function(x) x > 0)]) 
  }
}
print(cn)
paste(cn, collapse = "+")
glm.out = glm(paste("LogACG~",paste(cn,collapse = "+"),sep = ""), family = gaussian, data = salary)
cv.err = cv.glm(salary, glm.out, K = 10)$delta[1]

我的问题在于 for 循环。我试图构建一个循环，它将“vars”中的值附加到“model.use”中，但我似乎无法让它读取矩阵中的第二行。有什么建议么？谢谢

Answer 1

这里似乎发生了几件事。

您已将 LogACG 设置为倒数第二行 ("LogACG~，但它也是 vars 之一，由于 vars <- colnames(salary[c(2,3,7)])，它最终以 model.vars 结尾，所以这是不对的。

接下来你的第二个for 循环应该遍历model.spec.matrix 的行，即

for(i in 1:nrow(model.spec.matrix)){

并以编程方式捕获您可以执行的该行指示的列名（变量）

cn <- colnames(model.spec.matrix[sapply(model.spec.matrix[i,], function(x) x > 0)])

在循环中。您还应该将glm 和cv.glm 移动到循环内。

但是，这将覆盖 glm.out 和 cv.err 每次，因此您需要将它们创建为空列表并在每次迭代中附加列表。

所以最终产品将如下所示：

# Since you can't use LogACG to explain itself, 
#    suppose you meant to use Engineering as a candidate X
vars <- colnames(salary[c(2,3,8)])

# Make your grid
model.spec.matrix        <- expand.grid(list.to.expand)
names(model.spec.matrix) <- vars

glm.out <- list(rep(NA, nrow(model.spec.matrix)))
cv.err  <- list(rep(NA, nrow(model.spec.matrix)))

for(i in 1:nrow(model.spec.matrix)){
  cn <- colnames(model.spec.matrix[sapply(model.spec.matrix[i,], function(x) x > 0)])
  if(length(cn) == 0){cn <- "."}  
  tmp        <- glm(as.formula(paste("LogACG~",paste(cn,collapse = "+"),sep = "")), family = gaussian, data = salary)
  glm.out[i] <- capture.output(tmp$formula)
}
# > glm.out
# [[1]]
# [1] "LogACG ~ ."
# 
# [[2]]
# [1] "LogACG ~ Rank_Code"
# 
# [[3]]
# [1] "LogACG ~ Gender"
# 
# [[4]]
# [1] "LogACG ~ Rank_Code + Gender"
# 
# [[5]]
# [1] "LogACG ~ Engineering"
# 
# [[6]]
# [1] "LogACG ~ Rank_Code + Engineering"
# 
# [[7]]
# [1] "LogACG ~ Gender + Engineering"
# 
# [[8]]
# [1] "LogACG ~ Rank_Code + Gender + Engineering"

要获取列表元素中的整个模型对象，请替换

glm.out[i] <- capture.output(tmp$formula)

与

glm.out    <- append(glm.out, tmp)

即

for(i in 1:nrow(model.spec.matrix)){
  cn <- colnames(model.spec.matrix[sapply(model.spec.matrix[i,], function(x) x > 0)])
  if(length(cn) == 0){cn <- "."}  
  tmp        <- glm(as.formula(paste("LogACG~",paste(cn,collapse = "+"),sep = "")), family = gaussian, data = salary)
  tmp1       <- cv.glm(salary, tmp, K = 10)$delta[1]
  glm.out    <- append(glm.out, tmp)
  cv.err     <- append(cv.err, tmp1)
}

> tail(cv.err)
[[1]]
[1] 2.751025

[[2]]
[1] 2.758954

[[3]]
[1] 2.735063

[[4]]
[1] 2.768075

[[5]]
[1] 2.774433

[[6]]
[1] 2.748291

此外，您最好使用bestglm 之类的包，或仅使用默认包stats 中的step 函数（请参阅?step）。