我怎样才能组合几个变量并根据回归保持特定的集合

Question

我尝试制作如下示例数据：

set.seed(1)            # for reproducible example
x <- sample(100*20)
x <- matrix(x, nc = 20)     # 20 predictor
y <- 1 + 2*x[,1] + 3*x[,2] + 4*x[,3] + 5*x[,7] + 6*x[,8] + 7*x[,9] + rnorm(100)  # y depends on variables 1,2,3,7,8,9 only

df <- data.frame(y, as.matrix(x))

现在，我想组合 4 列 x 并保留 lm 模型的 R 高于 0.8 的所有这些组合

要在 Y 和 X 的 4 个变量之间建立模型，例如可以使用

fit = lm(Y~.,data=df[,c(2:6)])

我希望这 20 列中 4 个变量的所有组合的 R 回归高于 0.8

有人可以评论吗？

Answer 1

根据我的评论，我建议使用leaps 包，该包提供了一种算法，可以详尽地测试模型公式的每个变量组合，并返回一些指标（R 平方、BIC 等）。

您可以处理结果以获取符合您标准的变量列表（这里我采用 0.85 限制来获得较小的模型列表）。首先，拟合模型并指定 4 个变量的限制，我们将只保留 20 个最好的模型（总共有 19380 个可能的 4 个变量模型......）：

library(leaps)
fit <- regsubsets(y~., df, nvmax=4, nbest=20)

根据 r 平方限制（存储在摘要的另一个输出中）对输出表进行子集（它是一个布尔表，对于模型中保存的每个变量都具有 TRUE/FALSE）：

mytable <- data.frame(tail(summary(fit)$which, 20)[which(tail(summary(fit)$rsq, 20)>0.85), ])

排列它以获得最佳模型的变量名，以转置格式，第一个是 R² 中最好的：

output <- t(apply(mytable, 1, function(x) names(mytable)[x]))
####      [,1]           [,2] [,3] [,4] [,5] 
#### [1,] "X.Intercept." "X3" "X7" "X8" "X9" 
#### [2,] "X.Intercept." "X2" "X7" "X8" "X9" 
#### [3,] "X.Intercept." "X1" "X7" "X8" "X9" 
#### [4,] "X.Intercept." "X7" "X8" "X9" "X15"
#### [5,] "X.Intercept." "X4" "X7" "X8" "X9" 

如果您需要使用其中一种模型进行拟合，您可以像这样检索公式：

as.formula(paste("y ~ ", paste(output[1, -1], collapse=" + ")))
#### y ~ X3 + X7 + X8 + X9

或者干脆使用reformulate，感谢@Ben Bolker 的建议：

reformulate(output[1, -1], response="y")

编辑：使用 lasso 回归建模您的数据。

我使用的是 Hastie&Tishirani 改编的脚本，你还必须加载这些辅助函数here。我建议你首先掌握这项技术。首先，我创建数据并拆分训练集：

library(glmnet); set.seed(6)
train.ratio <- 0.75 
x      <- model.matrix(y~., df)[, -1]           
y      <- df$y
train.ind  <-   sample(1:nrow(x), floor(train.ratio * nrow(x)))
x.train    <-   x[ train.ind, ]; y.train    <-   y[ train.ind  ]
x.test     <-   x[-train.ind, ]; y.test     <-   y[-train.ind  ]
n          <-   nrow(x); n.train    <-   nrow(x.train); n.test     <-   nrow(x.test)

然后我进行模型校准并计算测试误差

grid <- 10^seq(4, -2, length=100) # increase range if needed

lasso.mod <- glmnet(x.train, y.train, alpha=1, lambda=grid) 
err.lasso <- 1/n.test * colSums((y.test - lasso.mod$a0[1] - x.test %*% lasso.mod$beta)^2)

绘制不同的结果

par(mfrow = c(2, 2))
frac.lasso <- plot.path(t(lasso.mod$beta), err = err.lasso)
plot.coef(t(lasso.mod$beta), lasso.mod$lambda, err.lasso)
plot.err(err.lasso, frac.lasso)
plot.err(err.lasso, lasso.mod$lambda)
par(mfrow = c(1, 1))

最终找出最好的模型是什么（很酷，它和以前一样！但不同的系数）。

best <- lasso.mod$beta[, which.min(err.lasso)]; best[best!=0]
####        X3        X7        X8        X9 
#### 0.2572322 1.4933962 1.7868181 2.4447500 

您也可以询问所有包含 4 个变量的案例（这里它们是相同的，只是 lambda 值和 coefs 不同）

lasso.mod$beta[, which(colSums(lasso.mod$beta!=0)==4)]