r - 在 R 中模拟数据集以进行模型选择

Question

我编写了一个代码来模拟 R 中的数据集，以了解反向选择在机器学习中的工作原理。我生成了 poly() 函数来编写多项式函数，然后想使用 Cp、BIC、调整后的 R^2 选择合适的多项式。

代码是:

 ###Generating dataset
 set.seed(1)
 X = rnorm(100)
 eps = rnorm(100)

 beta0 = 3
 beta1 = 2
 beta2 = -3
 beta3 = 0.3
 Y = beta0 + beta1 * X + beta2 * X^2 + beta3 * X^3 + eps
 
library(leaps)
data.full = data.frame(y = Y, x = X)
mod.full = regsubsets(y ~ poly(x, 10, raw = T), data = data.full, nvmax = 10)
mod.summary = summary(mod.full)

### Find the model size for best cp, BIC and adjr2
which.min(mod.summary$cp)

对于 cp、BIC 和调整后的 R^2，我得到的模型应该是多项式 3

但是，现在我想模拟 100 个数据集，看看在多少个数据集中我得到了正确的模型。我模拟了 100 个数据集，但现在我没有得到每个度量的多项式 3。我不太明白我做错了什么。我的模拟代码是:

###Generating 100 datasets
data <- replicate(100, rnorm(n=100))
epsilon <- replicate(100,rnorm(n=100))

###Formula (same as before)
 Y = beta0 + beta1 * data + beta2 * data^2 + beta3 * data^3 + epsilon
 data.full = data.frame(y = Y, x = data)

 ###Using polynomial terms
 mod.bwd = regsubsets(data.full$y.1 ~ poly(data.full$x.1, 10, raw = T), data = data.full, nvmax = 10, 
method = "backward")
bwd.summary = summary(mod.bwd)

which.min(bwd.summary$cp)
which.min(bwd.summary$bic)
which.max(bwd.summary$adjr2)

对于给定的子集 cp、Bic、adjr2 给出了不同的结果。例如，使用 y.1 和 x.1(模拟中的第一个数据集)给出以下结果:

which.min(bwd.summary$cp): 7

which.min(bwd.summary$bic): 4

which.max(bwd.summary$adjr2): 9

有人可以帮助我模拟这 100 个数据集时做错了什么吗？

Answer 1

如果我正确地阅读了您的代码，您可以在同一个模拟数据集上运行模型 100 次，而不是在所有 100 个模拟数据集上运行模型，这应该可以解决问题:

set.seed(42)

###Generating 100 datasets
data <- replicate(100, rnorm(n=100))
epsilon <- replicate(100,rnorm(n=100))

###Formula (same as before)
 Y = beta0 + beta1 * data + beta2 * data^2 + beta3 * data^3 + epsilon
 data.full = data.frame(y = Y, x = data)


res <- lapply(1:100, function(i){
 ###Using polynomial terms
 mod.bwd = regsubsets(data.full[[i]] ~ poly(data.full[[100+i]], 10, raw = T), data = data.full, nvmax = 10, 
method = "backward")
bwd.summary = summary(mod.bwd)

c(
which.min(bwd.summary$cp),
which.min(bwd.summary$bic),
which.max(bwd.summary$adjr2)
)
})

res <- do.call(rbind, res)

使用此 rng-seed，这会给出一些行，其中所有 cirteria 选择正确的模型。