Запуск регрессионной модели более 30, заданных c set.seed автоматически с использованием R - PullRequest
Новая
       21

Запуск регрессионной модели более 30, заданных c set.seed автоматически с использованием R

5 голосов
/ 28 марта 2020

Я работаю с несколькими моделями линейной регрессии.

Я хочу запустить модель линейной регрессии с разными 30 set.seed

Для пояснения, я делюсь кодом только с двумя моделями регрессии и 10 set.seed (В моем проекте у меня 12 регрессионных моделей, и каждая должна работать с 30 различными set.seeds)

Мне нужно решение, которое я могу запустить 30 set.seed для одной линейной регрессии модель, так что я могу go от моего ноутбука в течение периода работы (30 set.seeds). Затем я проделал то же самое для второй модели регрессии.

Есть ли способ автоматически выполнить код для 30 различных set.seed. Таким образом, я получил результат для каждого set.seed.

Я надеюсь, что все ясно, и я рад разъяснить больше.

NOTE

Bear Имейте в виду, что у меня есть четыре связанных блока с каждой моделью регрессии. Поэтому любое изменение с помощью set.seed или creatFolds может повлиять на другие блоки.

EDIT1

Используемый набор данных используется 

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)

1 Ответ

2 голосов
/ 02 апреля 2020

Если я вас правильно понимаю, вы хотите регрессировать V6 по всем другим переменным, используя как OLS, так и модель LAD. Вы хотите выбрать k=30 случайные "сгибы" с помощью createFolds и повторить процесс также n=30 раз. В результате вы хотите, чтобы отклонения для каждого повторения и каждого коэффициента.

Я бы обернул подходящую часть в функцию FX. Создайте n=30 семян с sample, l oop над ним с 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)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...