#r #for-loop #dplyr
#r #цикл for #dplyr
Вопрос:
Я пытаюсь смоделировать подобранные значения для изменения моделей разных подгрупп моих данных, которые снова основаны на случайной выборке другого подмножества моего исходного фрейма данных (минимальный пример, который я составил для этого вопроса, игнорирует случайную выборку и т.д., Что приводит к идентичным подобранным значениям для всех симуляций, но это не имеет значения). Я написал dplyr-код для хранения модели для каждой группы, создания новых значений x для прогнозирования установленных значений, их прогнозирования и т.д. и т.п. В результате получается один столбец с подобранными значениями, именно такими, как я хочу. Тем не менее, я хотел бы пройти весь процесс 1000 раз. Я, конечно, могу сделать это с помощью цикла for (как сделано ниже в примере), но есть ли возможность сделать такую вещь в dplyr-pipe line? И, может быть, немного ускорить весь процесс (мой исходный набор данных довольно большой, и цикл for занимает вечность)?
# making up data
dat <- data.frame("species" = seq(1:20), "col_A" = runif(20, min=1000, max=2500), "col_B" = runif(20, min = 0, max = 1500),
"maximum" = rep(2500, 20), "minimum" = rep(1000, 20), groups = rep(LETTERS[1:5], each = 4))
# functions to use with purrr
linear_mod <- function(dat) {
lm(col_A ~ col_B, data = dat)
}
# define parametres and an empty data frame to use in the for loop
runs <- 10
fitted_sim <- data.frame(matrix(data=NA,nrow=20,ncol=runs 1,byrow=FALSE)) #empty dataframe to contain fitted values for each alt
names(fitted_sim) <- as.factor(seq(1:runs 1))
# for-loop around my dplyr-code
for (j in 1:runs){
simul <- dat %>%
group_by(groups) %>%
nest(data = c(col_A, col_B, species)) %>%
mutate(model = map(data, linear_mod), # add model for every group
sim_data = list(seq(minimum, maximum, by = 10))) %>% # define new x-values for later predictions
unnest(sim_data) %>%
nest(sim_data = sim_data) %>%
mutate(fitted = map2(model, sim_data, ~predict(.x, col_B = sim_data, type = "response")), # predict values
unnest(fitted) # unnest predicted values to save in data frame
# save newly fitted values in fitted_sim data frame
fitted_sim[,1] <- simul$groups
fitted_sim[,1 j] <- simul$fitted
}
Спасибо за каждую подсказку!
РЕДАКТИРОВАТЬ: Это расширенный пример кода, включающий случайную выборку, описанную выше, но опущенную в моем первом примере:
# for-loop around my dplyr-code
for (j in 1:runs) {
simul <- dat %>%
group_by(groups) %>%
rowwise() %>%
mutate(vector_column = case_when(abs(col_A) == (maximum-minimum) ~ list(col_B),
sign(col_A) == 1 ~ list(dat$col_B[dat$col_B <= maximum - col_A]), # using list function to store vectors in a data.frame
sign(col_A) != 1 ~ list(dat$col_B[dat$col_B >= minimum col_A])),
helper = !is_empty(vector_column), # in case some of the vectors are empty so it is not possible to use sample
col_B_new = ifelse(helper, sample(vector_column, 1), NA),
helper = NULL,
sim_data = list(seq(minimum, maximum, by = 10))) %>% # define new x-values for later predictions
ungroup() %>% # get rid of rowwise()
group_by(groups) %>%
unnest(sim_data) %>%
nest(sim_data = sim_data) %>%
nest(data = c(col_A, col_B, col_B_new, species)) %>%
mutate(model = map(data, linear_mod),
fitted = map2(model, sim_data, ~predict(.x, col_B_new = sim_data, type = "response"))) %>% # predict values
unnest(fitted) # unnest predicted values to save in data frame
# save newly fitted values in fitted_sim data frame
fitted_sim[,1] <- simul$groups
fitted_sim[,1 j] <- simul$fitted
}
Ответ №1:
Не уверен, что именно вам понадобится из ваших simul data.frame в будущем, но отдельные элементы цикла кажутся независимыми друг от друга, поэтому вы могли бы использовать parallel:mclappy для их параллельного запуска. В моем примере я просто использую lapply , но это уже немного быстрее — если это вам как-то поможет…
library(tidyverse)
library(data.table)
library(microbenchmark)
# making up data
dat <- data.frame("species" = seq(1:20), "col_A" = runif(20, min=1000, max=2500), "col_B" = runif(20, min = 0, max = 1500),
"maximum" = rep(2500, 20), "minimum" = rep(1000, 20), groups = rep(LETTERS[1:5], each = 4))
# functions to use with purrr
linear_mod <- function(dat) {
lm(col_A ~ col_B, data = dat)
}
# define parametres and an empty data frame to use in the for loop
runs <- 10
fitted_sim <- data.frame(matrix(data=NA,nrow=20,ncol=runs 1,byrow=FALSE)) #empty dataframe to contain fitted values for each alt
names(fitted_sim) <- as.factor(seq(1:runs 1))
# I just wrapped your code in a function for benchmarking
run.simul <- function(){
for (j in 1:runs){
simul <- dat %>%
group_by(groups) %>%
nest(data = c(col_A, col_B, species)) %>%
mutate(model = map(data, linear_mod), # add model for every group
sim_data = list(seq(minimum, maximum, by = 10))) %>% # define new x-values for later predictions
unnest(sim_data) %>%
nest(sim_data = sim_data) %>%
mutate(fitted = map2(model, sim_data, ~predict(.x, col_B = sim_data, type = "response"))) %>% # predict values
unnest(fitted) # unnest predicted values to save in data frame
# save newly fitted values in fitted_sim data frame
fitted_sim[,1] <- simul$groups
fitted_sim[,1 j] <- simul$fitted
}
fitted_sim
}
# my version (sorry for using `data.table`...)
Dat <- data.table(dat, key="groups")
getFits <- function(x){
x[, fitted := predict(lm(col_A ~ col_B), .(list(seq(minimum[1], maximum[1], by = 10)))), by = groups]
x[, .(species, groups, fitted)]
}
# get the results into a data.frame; use mclapply instead of lapply if you like
dcast(rbindlist(lapply(seq_len(runs), function(z) getFits(data.table(dat, key="groups"))), idcol = "run"), ... ~ run, value.var = "fitted")
#> species groups 1 2 3 4 5 6
#> 1: 1 A 1767.621 1767.621 1767.621 1767.621 1767.621 1767.621
#> 2: 2 A 1376.679 1376.679 1376.679 1376.679 1376.679 1376.679
#> 3: 3 A 1523.100 1523.100 1523.100 1523.100 1523.100 1523.100
#> 4: 4 A 1794.593 1794.593 1794.593 1794.593 1794.593 1794.593
#> 5: 5 B 1272.408 1272.408 1272.408 1272.408 1272.408 1272.408
#> 6: 6 B 1967.709 1967.709 1967.709 1967.709 1967.709 1967.709
#> 7: 7 B 1792.934 1792.934 1792.934 1792.934 1792.934 1792.934
#> 8: 8 B 2318.699 2318.699 2318.699 2318.699 2318.699 2318.699
#> 9: 9 C 2017.187 2017.187 2017.187 2017.187 2017.187 2017.187
#> 10: 10 C 1899.827 1899.827 1899.827 1899.827 1899.827 1899.827
#> 11: 11 C 1734.953 1734.953 1734.953 1734.953 1734.953 1734.953
#> 12: 12 C 1834.046 1834.046 1834.046 1834.046 1834.046 1834.046
#> 13: 13 D 1797.615 1797.615 1797.615 1797.615 1797.615 1797.615
#> 14: 14 D 1915.832 1915.832 1915.832 1915.832 1915.832 1915.832
#> 15: 15 D 1841.489 1841.489 1841.489 1841.489 1841.489 1841.489
#> 16: 16 D 1798.442 1798.442 1798.442 1798.442 1798.442 1798.442
#> 17: 17 E 1641.359 1641.359 1641.359 1641.359 1641.359 1641.359
#> 18: 18 E 1631.253 1631.253 1631.253 1631.253 1631.253 1631.253
#> 19: 19 E 1634.197 1634.197 1634.197 1634.197 1634.197 1634.197
#> 20: 20 E 1634.991 1634.991 1634.991 1634.991 1634.991 1634.991
#> 7 8 9 10
#> 1: 1767.621 1767.621 1767.621 1767.621
#> 2: 1376.679 1376.679 1376.679 1376.679
#> 3: 1523.100 1523.100 1523.100 1523.100
#> 4: 1794.593 1794.593 1794.593 1794.593
#> 5: 1272.408 1272.408 1272.408 1272.408
#> 6: 1967.709 1967.709 1967.709 1967.709
#> 7: 1792.934 1792.934 1792.934 1792.934
#> 8: 2318.699 2318.699 2318.699 2318.699
#> 9: 2017.187 2017.187 2017.187 2017.187
#> 10: 1899.827 1899.827 1899.827 1899.827
#> 11: 1734.953 1734.953 1734.953 1734.953
#> 12: 1834.046 1834.046 1834.046 1834.046
#> 13: 1797.615 1797.615 1797.615 1797.615
#> 14: 1915.832 1915.832 1915.832 1915.832
#> 15: 1841.489 1841.489 1841.489 1841.489
#> 16: 1798.442 1798.442 1798.442 1798.442
#> 17: 1641.359 1641.359 1641.359 1641.359
#> 18: 1631.253 1631.253 1631.253 1631.253
#> 19: 1634.197 1634.197 1634.197 1634.197
#> 20: 1634.991 1634.991 1634.991 1634.991
run.simul()
#> 1 2 3 4 5 6 7 8 9
#> 1 A 1767.621 1767.621 1767.621 1767.621 1767.621 1767.621 1767.621 1767.621
#> 2 A 1376.679 1376.679 1376.679 1376.679 1376.679 1376.679 1376.679 1376.679
#> 3 A 1523.100 1523.100 1523.100 1523.100 1523.100 1523.100 1523.100 1523.100
#> 4 A 1794.593 1794.593 1794.593 1794.593 1794.593 1794.593 1794.593 1794.593
#> 5 B 1272.408 1272.408 1272.408 1272.408 1272.408 1272.408 1272.408 1272.408
#> 6 B 1967.709 1967.709 1967.709 1967.709 1967.709 1967.709 1967.709 1967.709
#> 7 B 1792.934 1792.934 1792.934 1792.934 1792.934 1792.934 1792.934 1792.934
#> 8 B 2318.699 2318.699 2318.699 2318.699 2318.699 2318.699 2318.699 2318.699
#> 9 C 2017.187 2017.187 2017.187 2017.187 2017.187 2017.187 2017.187 2017.187
#> 10 C 1899.827 1899.827 1899.827 1899.827 1899.827 1899.827 1899.827 1899.827
#> 11 C 1734.953 1734.953 1734.953 1734.953 1734.953 1734.953 1734.953 1734.953
#> 12 C 1834.046 1834.046 1834.046 1834.046 1834.046 1834.046 1834.046 1834.046
#> 13 D 1797.615 1797.615 1797.615 1797.615 1797.615 1797.615 1797.615 1797.615
#> 14 D 1915.832 1915.832 1915.832 1915.832 1915.832 1915.832 1915.832 1915.832
#> 15 D 1841.489 1841.489 1841.489 1841.489 1841.489 1841.489 1841.489 1841.489
#> 16 D 1798.442 1798.442 1798.442 1798.442 1798.442 1798.442 1798.442 1798.442
#> 17 E 1641.359 1641.359 1641.359 1641.359 1641.359 1641.359 1641.359 1641.359
#> 18 E 1631.253 1631.253 1631.253 1631.253 1631.253 1631.253 1631.253 1631.253
#> 19 E 1634.197 1634.197 1634.197 1634.197 1634.197 1634.197 1634.197 1634.197
#> 20 E 1634.991 1634.991 1634.991 1634.991 1634.991 1634.991 1634.991 1634.991
#> 10 NA
#> 1 1767.621 1767.621
#> 2 1376.679 1376.679
#> 3 1523.100 1523.100
#> 4 1794.593 1794.593
#> 5 1272.408 1272.408
#> 6 1967.709 1967.709
#> 7 1792.934 1792.934
#> 8 2318.699 2318.699
#> 9 2017.187 2017.187
#> 10 1899.827 1899.827
#> 11 1734.953 1734.953
#> 12 1834.046 1834.046
#> 13 1797.615 1797.615
#> 14 1915.832 1915.832
#> 15 1841.489 1841.489
#> 16 1798.442 1798.442
#> 17 1641.359 1641.359
#> 18 1631.253 1631.253
#> 19 1634.197 1634.197
#> 20 1634.991 1634.991
# benchmarking
microbenchmark(
a=dcast(rbindlist(lapply(seq_len(runs), function(z) getFits(data.table(dat, key="groups"))), idcol = "run"), ... ~ run, value.var = "fitted"),
b=run.simul(), times= 10L
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval cld
#> a 46.46851 48.65027 51.09618 49.88775 54.7177 56.2888 10 a
#> b 389.76893 410.57588 418.16993 415.74493 424.2042 448.4604 10 b
Создано 2020-08-13 пакетом reprex (версия 0.3.0)
Комментарии:
1. Это здорово, спасибо! Я не очень хорошо знаком с таблицей данных, но я определенно должен в нее заглянуть. Однако для моих реальных данных для каждого раунда цикла for я случайным образом выбираю одно новое значение col_B для каждого вида из столбца списка в моем фрейме данных, может ли это быть интегрировано в ваш код? Работает ли таблица данных со столбцами списка?
2. Я добавил расширенную версию примера скрипта, чтобы проиллюстрировать, что я имею в виду. Не стесняйтесь игнорировать это, хотя вы прекрасно ответили на мой первый вопрос.
3. Если вы не изменили свою
linear_modфункцию, вы не изменили модели при редактировании. Вы хотите сделатьlm(col_A ~ col_B_new)? Кроме того, ваша выборка, похоже, выполняется по группам дляabs(col_A) == (maximum-minimum), но в противном случае по всем группам, это предназначено? И да,data.tableимеет дело со столбцами списка, как показано в моейgetFitsфункции, см. Строку с.(list(seq(minimum[1], maximum[1], by = 10)), где.()список переносится в другой список, чтобы сохранить внутренний список в столбце.4. Нет, я не менял модели — я просто не был уверен, как интегрировать дополнительные условия, но теперь я понял это. Большое спасибо за вашу помощь!