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

Simulating datasets in R for model selection

我编写了一个代码来模拟 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 个数据集时我做错了什么。

如果我没看错你的代码,你 运行 同一个模拟数据集上的模型 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 select 都是正确的模型。