我正在使用几个线性回归模型。
我想用不同的 30 set.seed 运行线性回归模型
为了澄清起见，我只与两个回归模型和 10 个 set.seed 共享代码(在我的项目中，我有 12 个回归模型，每个模型都应该使用 30 个不同的 set.seeds 运行)
我需要一个解决方案，我可以为一个线性回归模型运行 30 set.seed，这样我就可以在运行期间离开我的笔记本电脑(30 set.seeds)。然后我对第二个回归模型做了同样的事情。
有没有办法在 30 个不同的 set.seed 上自动运行代码。所以我得到了每个 set.seed 的结果。
我希望一切都清楚，我很乐意澄清更多。
注意
请记住，每个回归模型都有四个相关的块。因此，对 set.seed 或 creatFolds 的任何更改都可能影响其他块。
编辑1
dataset用过的
wdbc <- read.delim("airfoil_self_noise.dat",header=F)
wdbcc=as.data.frame(scale(wdbc))
#set.seed(1)
#set.seed(2)
#set.seed(3)
#set.seed(4)
...
k = 30
folds <- createFolds(wdbcc$V6, k = k, list = TRUE, returnTrain = TRUE)
## Ordinary Least Square regression ##
#Block A
lm = list()
for (i in 1:k) {
lm[[i]] = lm(V6~ ., data = wdbcc[folds[[i]],])
}
#Block B
lm_coef = list()
lm_coef_var = list()
for(j in 1:(lm[[1]]$coefficients %>% length())){
for(i in 1:k){
lm_coef[[i]] = lm[[i]]$coefficients[j]
lm_coef_var[[j]] = lm_coef %>% unlist() %>% var()
}
}
#Block C
lm_var = unlist(lm_coef_var)
lm_df = cbind(coefficients = lm[[1]]$coefficients%>% names() %>% as.data.frame()
, variance = lm_var %>% as.data.frame())
colnames(lm_df) = c("coefficients", "variance_lm")
#Block D
lm_var_sum = sum(lm_var)
PQSQ-回归
X=list()
Y=list()
for (i in 1:k) {
n=wdbcc[folds[[i]],-6]
m=wdbcc[folds[[i]],6]
X[[i]]=n
Y[[i]]=m
}
#Block A
lmPQSQ1 = list()
for (i in 1:k) {
lmPQSQ1[[i]] = PQSQRegression(X[[i]],Y[[i]],0.01,data = wdbcc[folds[[i]],])
}
lmmPQSQ1=list()
for (i in 1:k) {
L=list(coefficients = c(lmPQSQ1[[i]][[2]],lmPQSQ1[[i]][[1]]))
lmmPQSQ1[[i]]=L
}
#Block B
lm_coefPQSQ1 = list()
lm_coef_varPQSQ1 = list()
for(j in 1:(lmmPQSQ1[[1]]$coefficients %>% length())){
for(i in 1:k){
lm_coefPQSQ1[[i]] = lmmPQSQ1[[i]]$coefficients[j]
lm_coef_varPQSQ1[[j]] = lm_coefPQSQ1 %>% unlist() %>% var()
}
}
#Block C
lm_varPQSQ1 = unlist(lm_coef_varPQSQ1)
lm_dfPQSQ1 = variance = lm_varPQSQ1 %>% as.data.frame()
#Block D
PQSQ1_var_sum = sum(lm_varPQSQ1)
最佳答案
如果我理解正确，你想回归 V6使用 OLS 和 LAD 模型的所有其他变量。您要选择 k=30使用 createFolds 的随机“折叠”并重复此过程 n=30次。因此，您需要每次重复和每个系数的方差。
我会将拟合部分包装成一个函数 FX .生成 n=30种子与 sample , 用 lapply 循环重复 FX n=30次。
FX <- function(seed, data, k=30) {
set.seed(seed) ## sets seed for each iteration
folds <- createFolds(data[, "V6"], k=k, list=TRUE, returnTrain=TRUE) ## folds
## OLS
lm1 <- lapply(folds, function(folds) lm(V6 ~ ., data=data[folds, ]))
lm.coefs <- t(sapply(lm1, coef)) ## lm coefficients
## LAD
lad1 <- lapply(folds, function(folds) lad(V6 ~ ., data=data[folds, ], method="BR"))
lad.coefs <- t(sapply(lad1, coef)) ## lad coefficients
## calculate column variances for both coef matrices
## use `attr<-` to add the seed as an attribute if you want
return(`attr<-`(cbind(lm=apply(lm.coefs, 2, var), lad=apply(lad.coefs, 2, var)),
"seed", seed))
}
seeds <- 1:30 ## specific seeds 1, 2, ... 30
## if you want non-consecutive specific seeds, do:
# set.seed(42) ## set some initial seed
# n <- 30 ## n. o. seeds
# seeds <- sample(1:1e6, n) ## sample seeds for `FX`
res <- lapply(seeds, FX, data=wdbcc) ## lapply loop
结果
这会生成一个长度为 30 的列表，其中包含每个重复、每个模型和每个系数的方差矩阵。
res[1:2] ## first two lists
# [[1]]
# lm lad
# (Intercept) 9.104280e-06 1.273920e-05
# V1 2.609623e-05 6.992753e-05
# V2 7.082099e-05 2.075875e-05
# V3 1.352299e-05 1.209651e-05
# V4 7.986000e-06 9.273005e-06
# V5 5.545298e-05 1.535849e-05
# attr(,"seed")
# [1] 1
#
# [[2]]
# lm lad
# (Intercept) 4.558722e-06 2.031707e-05
# V1 2.256583e-05 9.291900e-05
# V2 6.519648e-05 2.768443e-05
# V3 1.800889e-05 9.983524e-06
# V4 1.131813e-05 1.174496e-05
# V5 3.866105e-05 1.022452e-05
# attr(,"seed")
# [1] 2
length(res)
# [1] 30
要计算每个种子的方差总和，您可以使用
colSums在 sapply .# sum of variances
sov <- t(sapply(res, colSums))
dim(sov)
# [1] 30 2
head(sov)
# lm lad
# [1,] 1.829835e-04 0.0001401535
# [2,] 1.603091e-04 0.0001728735
# [3,] 1.003093e-04 0.0001972869
# [4,] 1.460591e-04 0.0001508251
# [5,] 9.915082e-05 0.0001262106
# [6,] 1.425996e-04 0.0001478449
了解
lapply 的迭代是什么确实，考虑一下:## provide the values of first iteration for arguments of function `FX`
seed <- 1
data <- wdbcc
k <- 30
## first iteration of `lapply`
set.seed(seed)
folds <- createFolds(data[, "V6"], k=k, list=TRUE, returnTrain=TRUE) ## folds
## OLS
lm1 <- lapply(folds, function(folds) lm(V6 ~ ., data=data[folds, ]))
lm.coefs <- t(sapply(lm1, coef)) ## lm coefficients
dim(lm.coefs)
# [1] 30 6
head(lm.coefs)
# (Intercept) V1 V2 V3 V4 V5
# Fold01 -0.0039130125 -0.5806272 -0.3564769 -0.4804492 0.2271908 -0.2805472
# Fold02 0.0013260444 -0.5863764 -0.3533327 -0.4759213 0.2253128 -0.2874691
# Fold03 0.0006791787 -0.5890787 -0.3678586 -0.4832066 0.2220979 -0.2739124
# Fold04 -0.0010721593 -0.5868079 -0.3722466 -0.4895328 0.2227811 -0.2749657
# Fold05 0.0021856620 -0.5850165 -0.3495360 -0.4810657 0.2235410 -0.2936287
# Fold06 0.0001486909 -0.5872607 -0.3677774 -0.4848523 0.2275780 -0.2823764
## LAD (same as OLS)
lad1 <- lapply(folds, function(folds) lad(V6 ~ ., data=data[folds, ], method="BR"))
lad.coefs <- t(sapply(lad1, coef)) ## lad coefficients
## return, throws variances for each coefficient of each model in a matrix
## the seed is added as an attribute, to be able to identify it later
res.1 <- `attr<-`(cbind(var.lm=apply(lm.coefs, 2, var),
var.lad=apply(lad.coefs, 2, var)),
"seed", seed)
res.1
# var.lm var.lad
# (Intercept) 9.104280e-06 1.273920e-05
# V1 2.609623e-05 6.992753e-05
# V2 7.082099e-05 2.075875e-05
# V3 1.352299e-05 1.209651e-05
# V4 7.986000e-06 9.273005e-06
# V5 5.545298e-05 1.535849e-05
# attr(,"seed")
# [1] 1
比较
res.1与列表的第一个元素 res以上。sov.1 <- colSums(res.1)
sov.1
# var.lm var.lad
# 0.0001829835 0.0001401535
比较
sov.1与矩阵的第一行 sov以上。编辑
对于 带矩阵表示法的回归函数 ，例如
lm.fit ，我们可以使用 model.matrix并预先进行子集化，见行 lm2.coefs在函数中；比较 lm和 lm2 res2 中的列下面，他们是平等的。 ( lm.fit 也比 lm 快，因为它省略了不必要的计算，你只需要系数；因此你实际上可以用 lm 行替换 lm.fit 。也可能有一种方法与 lad 使用lsfit 在代码中，但老实说，我对 lad 太陌生，无法为您提供此解决方案。)另请注意，为了简洁起见，我使用
sapply 将每个模型的两条线合并为一条线。直接上$coefficients . sapply作品为 lapply但抛出一个矩阵；请注意，我们需要t逃跑。FX2 <- function(seed, data, k=30) {
set.seed(seed) ## sets seed for each iteration
folds <- createFolds(data[, "V6"], k=k, list=TRUE, returnTrain=TRUE) ## draw folds
lm.coefs <- t(sapply(folds, function(f) lm(V6 ~ ., data=data[f, ])$coef))
lm2.coefs <- t(sapply(folds, function(f) {
data2 <- data[f, ]
lm.fit(x=model.matrix(V6 ~ ., data2), y=data2[,"V6"])$coef
}))
lad.coefs <- t(sapply(folds, function(f) lad(V6 ~ ., data=data[f, ], method="BR")$coef))
return(`attr<-`(cbind(lm=apply(lm.coefs, 2, var),
lm2=apply(lm2.coefs, 2, var),
lad=apply(lad.coefs, 2, var)),
"seed", seed))
}
seeds <- 1:30
res.2 <- lapply(seeds, FX2, data=wdbcc) ## lapply loop
res.2[1:2]
# [[1]]
# lm lm2 lad
# (Intercept) 9.104280e-06 9.104280e-06 1.273920e-05
# V1 2.609623e-05 2.609623e-05 6.992753e-05
# V2 7.082099e-05 7.082099e-05 2.075875e-05
# V3 1.352299e-05 1.352299e-05 1.209651e-05
# V4 7.986000e-06 7.986000e-06 9.273005e-06
# V5 5.545298e-05 5.545298e-05 1.535849e-05
# attr(,"seed")
# [1] 1
#
# [[2]]
# lm lm2 lad
# (Intercept) 4.558722e-06 4.558722e-06 2.031707e-05
# V1 2.256583e-05 2.256583e-05 9.291900e-05
# V2 6.519648e-05 6.519648e-05 2.768443e-05
# V3 1.800889e-05 1.800889e-05 9.983524e-06
# V4 1.131813e-05 1.131813e-05 1.174496e-05
# V5 3.866105e-05 3.866105e-05 1.022452e-05
# attr(,"seed")
# [1] 2
数据和图书馆:
invisible(lapply(c("caret", "L1pack"), library, character.only=TRUE))
wdbcc <- read.delim("airfoil_self_noise.dat", header=F)
wdbcc[] <- lapply(wdbcc, scale)
关于使用 R 自动运行超过 30 个特定 set.seed 的回归模型，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/60895347/