使用生成的随机数据重复一个过程，并将结果保存在 R 中的数据框中

Question

我像这样在 R 中创建随机数据：

data<-matrix(rnorm(100*5,mean=0,sd=1), 100, 5) 
colnames(data) <- c("X1", "X2", "X3", "X4", "X5")
data <- as.data.frame(data)
a <- 5 
b <- 0.8
c <- 100

然后我想“玩”这些数据的相关性并执行以下操作

data[,2] <- a*data[,1] - b*rnorm(c)
data[,3] <- a*data[,1] + b*rnorm(c)
data[,4] <- a*data[,1] - b*rnorm(c)

之后我执行以下代码

data<-matrix(rnorm(100*5,mean=0,sd=1), 100, 5) 
colnames(data) <- c("X1", "X2", "X3", "X4", "X5")
data <- as.data.frame(data)
a <- 5 
b <- 0.8
c <- 100

data[,2] <- a*data[,1] - b*rnorm(c)
data[,3] <- a*data[,1] + b*rnorm(c)
data[,4] <- a*data[,1] - b*rnorm(c)

library(glmnet)
library(coefplot)

A <- as.matrix(data)
set.seed(1)
results <- lapply(seq_len(ncol(A)), function(i) {
  list(
    cvfit = cv.glmnet(A[, -i] , A[, i] , standardize = TRUE , type.measure = "mse" , nfolds = 10 , alpha = 1)
  )
})

lam <- as.data.frame(`names<-`(
  lapply(results, function(x) (x$cvfit$lambda.min)), 
  paste0("X", seq_along(results))
))
    
    sigma<- matrix(rnorm(1*5,mean=0,sd=1), 1, 5) 
colnames(sigma) <- c("X1", "X2", "X3", "X4", "X5")
as.vector(sigma)
sub1.sigma <- subset(sigma, select = sigma <= sum(lam))
sub2.sigma <- subset(sigma, select = sigma <= 2*sum(lam))
sub3.sigma <- subset(sigma, select = sigma <= 3*sum(lam))

这会产生一个称为 sigma 的向量 1x5 和 3 个向量 sub1.sigma、sub2.sigma、sub3.sigma 像下面这样

   > sigma
      X1      X2        X3 X4 X5        
     38.64019 624.4896  0  0  0 
> sub1.sigma
        X1  X3  X4  X5       
1 38.64019  0   0   0 
 
> sub2.sigma
        X1  X3 X4 X5          
1 38.64019  0  0  0 

> sub3.sigma
        X1  X3 X4 X5         
1 38.64019  0  0  0 

生成的数据是随机的，我通常使用set.seed() 来产生相同的结果。如果可以不修改主代码，我想运行我的代码 100 次（每次使用不同的数据）并将相应的结果保存在 4 个数据帧中sigmasub1.sigma，sub2.sigma , sub3.sigma 以便比较它们。有没有办法在 R 中实现这一点？

基于 cmets，我设法创建了以下内容，但似乎仍然没有给出预期的结果。首先代码[1:10] 显示 10 个向量，它们代表什么？西格玛？这些是每次运行的 sigma 吗？我怎样才能让它也计算 sub.sigma？

set.seed(2021)

code <- replicate(10,{
data<-matrix(rnorm(100*5,mean=0,sd=1), 100, 5) 
colnames(data) <- c("X1", "X2", "X3", "X4", "X5")
data <- as.data.frame(data)
a <- 5 
b <- 0.8
c <- 100

data[,2] <- a*data[,1] - b*rnorm(c)
data[,3] <- a*data[,1] + b*rnorm(c)
data[,4] <- a*data[,1] - b*rnorm(c)

library(glmnet)
library(coefplot)

A <- as.matrix(data)
set.seed(1)
results <- lapply(seq_len(ncol(A)), function(i) {
  list(
    cvfit = cv.glmnet(A[, -i] , A[, i] , standardize = TRUE , type.measure = "mse" , nfolds = 10 , alpha = 1)
  )
})

lam <- as.data.frame(`names<-`(
  lapply(results, function(x) (x$cvfit$lambda.min)), 
  paste0("X", seq_along(results))
))

sigma<- matrix(rnorm(1*5,mean=0,sd=1), 1, 5) 
colnames(sigma) <- c("X1", "X2", "X3", "X4", "X5")
as.vector(sigma)
sub1.sigma <- subset(sigma, select = sigma <= sum(lam))
sub2.sigma <- subset(sigma, select = sigma <= 2*sum(lam))
sub3.sigma <- subset(sigma, select = sigma <= 3*sum(lam))
}, simplify = FALSE)

code[1:10]
sigmas <- as.data.frame(do.call(rbind,lapply(code, sigma)))

Answer 1

出于各种原因，我喜欢保留模型，因此我将从模型运行列表开始。在这种情况下，replicate(n, ..., simplify=FALSE) 返回一个 list 我们需要的任何东西。 （记录一下，相当于lapply(seq_len(n), function(ign) ...)。）

（旁注：我没有安装glmnet，所以我会用一个简单/荒谬的lm 模型来模仿它。）

set.seed(2021)
models <- replicate(10, {
  zany_numbers <- runif(32) # nrow(mtcars)
  lm(zany_numbers ~ mpg + disp + cyl, data = mtcars)
}, simplify = FALSE)
models[1:2]
# [[1]]
# Call:
# lm(formula = zany_numbers ~ mpg + disp + cyl, data = mtcars)
# Coefficients:
# (Intercept)          mpg         disp          cyl  
#    8.85e-02     2.01e-02    -5.29e-05     2.41e-02  
# [[2]]
# Call:
# lm(formula = zany_numbers ~ mpg + disp + cyl, data = mtcars)
# Coefficients:
# (Intercept)          mpg         disp          cyl  
#    -0.52302      0.02485     -0.00122      0.13071  

在这里，我们可以制作您想要的任何类型的框架。

coefs <- as.data.frame(do.call(rbind, lapply(models, coef)))
coefs
#    (Intercept)       mpg      disp     cyl
# 1       0.0885  0.020114 -5.29e-05  0.0241
# 2      -0.5230  0.024847 -1.22e-03  0.1307
# 3       0.0856  0.014215  1.18e-03 -0.0225
# 4       0.4899  0.012013  1.08e-03 -0.0876
# 5      -0.6926  0.024653 -1.16e-03  0.1499
# 6       0.0862  0.010497 -5.02e-04  0.0389
# 7       0.8358 -0.008419 -6.64e-04 -0.0141
# 8       0.3679 -0.000198  1.44e-03 -0.0391
# 9       0.4360  0.011994 -9.59e-05 -0.0303
# 10      0.2276  0.003659 -1.14e-03  0.0651

（您可能需要清理那里的名称。）

如果您愿意，可以将do.call(rbind, ...) 替换为data.table::rbindlist(...) 或dplyr::bind_rows(...)。

从这个models，并使用相同的框架列表do.call(rbind, ...)后续，您可以生成伴随框架，例如

otherstats <- as.data.frame(do.call(rbind, lapply(models, function(mdl) summary(mdl)[c("r.squared", "adj.r.squared")])))
otherstats
#    r.squared adj.r.squared
# 1      0.104       0.00745
# 2      0.144        0.0523
# 3      0.044       -0.0584
# 4      0.202         0.117
# 5      0.149        0.0573
# 6     0.0639       -0.0364
# 7     0.0586       -0.0422
# 8      0.137        0.0446
# 9      0.241          0.16
# 10    0.0814        -0.017