对不同值的组合应用 foreach 循环答案

【问题标题】：Applying a foreach loop over combinations of different values对不同值的组合应用 foreach 循环
【发布时间】：2021-02-01 15:23:19
【问题描述】：

我有一些数据想要应用某种蒙特卡罗模拟。我有 n 和 k 的不同组合，我想在这些不同组合上运行 forloop。以下代码适用于 n 和 k 的单个值

代码：

formula <- as.formula(predictions ~ age + daysInHospital)
result <- glm(formula = formula, data = df, family = binomial)

Tdays = 100 # number of days
k = 6 # k new people coming to the hospital each day - randomly taken from the df
N = 89 # number of available beds each day
# objective: decide who stays in the hospital and who is released based on the probabilities
# policy: keep in the hospital the N patients with the highest prob of mortality
# each day, we have a cut-off and patients with at least this prob of mortality will stay in hospital for another day
# i.e. newcomers will have 0 days in hospital but everybody else will have + 1 day in hospital
MCtot = 5 # total number of monte carlos simulations to run

set.seed(0)
out <- data.frame()

for(indMC in 1:MCtot){
  print(paste("Current MC run: ", indMC, "out of: ", MCtot ))
  # initial patients in hospital for this MC loop
  patients = df %>% 
    sample_n(N)
  
  for(t in 1:Tdays){
    # take newcomers, with different prob and with daysInHospital=0
    newcomers = df %>% 
      sample_n(k) %>% 
      mutate(
        daysInHospital = 0 # since these are the newcommers they will have 0 days in hospital
      )
    
    # apply policy: 'keep the N patients with highest prob'
    # make a pool with newcomers and patients
    aux1 = patients
    aux1 = rbind(aux1, newcomers, fill = TRUE)
    # make predictions at period t,
    aux1$t_pred = predict(result, aux1, type = "response")
    # take the N-th highest
    cut_off <- aux1 %>% 
      arrange(desc(t_pred)) %>% 
      slice(N) %>% 
      pull(t_pred)
    # all patients with prob >= cut_off will stay in hospital one more day
    patients <- aux1 %>% 
      filter(t_pred >= cut_off)
    
    avg_prob <- patients %>% 
      summarise(
        avg_prob = mean(t_pred)
      ) %>% 
      pull(avg_prob)
    
    # remove auxiliar
    aux1 <- NULL
    # store information of current day
    aux1 <- data.frame(
      MCrun = indMC,
      period = t,
      cut_off = cut_off,
      avg_prob = avg_prob,
      n = N,
      k = k
    )
    out <- rbind(out, aux1)
    # update daysinHospital in patients and remove auxiliar
    
    # update daysinHospital in patients
    patients <- patients %>% 
      mutate(
        daysInHospital + 1
      )
  }
}

现在我想生成 n 和 k 的随机组合，并应用相同的 forloop 将结果存储在 out 中。

我尝试将 forloop 包装在 foreach 循环中，但我没有得到正确的输出，它返回并出错。

输出：

[1] "Current k:  4 Current n:  90"
[1] "Current MC run:  1 out of:  5"
[1] "Current MC run:  2 out of:  5"
[1] "Current MC run:  3 out of:  5"
[1] "Current MC run:  4 out of:  5"
[1] "Current MC run:  5 out of:  5"
[1] "Current k:  96 Current n:  21"
[1] "Current MC run:  1 out of:  5"
[1] "Current k:  57 Current n:  4"
[1] "Current MC run:  1 out of:  5"
[1] "Current MC run:  2 out of:  5"
[1] "Current MC run:  3 out of:  5"
[1] "Current MC run:  4 out of:  5"
[1] "Current MC run:  5 out of:  5"
Error in { : 
  task 2 failed - "`size` must be less or equal than 94 (size of data), set `replace` = TRUE to use sampling with replacement."

代码（中断）

formula <- as.formula(predictions ~ age + daysInHospital)
result <- glm(formula = formula, data = df, family = binomial)

Tdays = 100 # number of days

K = floor(runif(3, min = 0, max = 100))
N = floor(runif(3, min = 0, max = 100))

MCtot = 5 # total number of monte carlos simulations to run

set.seed(0)
out <- data.frame()



foreach(k = K, n = N) %do% {
  print(paste("Current k: ", k, "Current n: ", n))
  for(indMC in 1:MCtot){
    print(paste("Current MC run: ", indMC, "out of: ", MCtot ))
    # initial patients in hospital for this MC loop
    patients = df %>% 
      sample_n(n)                                 # changed n here (from upper to lowercase)
    
    for(t in 1:Tdays){
      # take newcomers, with different prob and with daysInHospital=0
      newcomers = df %>% 
        sample_n(k) %>% 
        mutate(
          daysInHospital = 0 # since these are the newcommers they will have 0 days in hospital
        )
      
      # apply policy: 'keep the N patients with highest prob'
      # make a pool with newcomers and patients
      aux1 = patients
      aux1 = rbind(aux1, newcomers, fill = TRUE)
      # make predictions at period t,
      aux1$t_pred = predict(result, aux1, type = "response")
      # take the N-th highest
      cut_off <- aux1 %>% 
        arrange(desc(t_pred)) %>% 
        slice(n) %>%                                  # changed n here (from upper to lowercase) 
        pull(t_pred)
      # all patients with prob >= cut_off will stay in hospital one more day
      patients <- aux1 %>% 
        filter(t_pred >= cut_off)
      
      avg_prob <- patients %>% 
        summarise(
          avg_prob = mean(t_pred)
        ) %>% 
        pull(avg_prob)
      
      # remove auxiliar
      aux1 <- NULL
      # store information of current day
      aux1 <- data.frame(
        MCrun = indMC,
        period = t,
        cut_off = cut_off,
        avg_prob = avg_prob,
        n = n,                                  # changed n here (from upper to lowercase)
        k = k
      )
      out <- rbind(out, aux1)
      # update daysinHospital in patients and remove auxiliar
      
      # update daysinHospital in patients
      patients <- patients %>% 
        mutate(
          daysInHospital + 1
        )
    }
  }
}

关于如何将 forloop 应用于 n 和 k 上的不同组合的任何提示。

此外，我尝试重新创建 forloop，但使用了函数 - 如果需要，我可以放置函数的代码 - 它不完整，因为我无法在每次迭代时更新 daysInHospital 变量 patients <- patients %>% mutate(daysInHospital + 1)

包：

library(dplyr)
library(ggplot2)
library(tidyverse)
library(foreach)

数据：

    df <- structure(list(predictions = c(0.456592172384262, 0.311251223087311, 
0.322826206684113, 0.320436120033264, 0.420515507459641, 0.311251223087311, 
0.340740621089935, 0.344267100095749, 0.33494707942009, 0.316163510084152, 
0.439167380332947, 0.45088067650795, 0.348440110683441, 0.348440110683441, 
0.364672362804413, 0.311251223087311, 0.311251223087311, 0.31931933760643, 
0.363355785608292, 0.311251223087311, 0.320436120033264, 0.419657677412033, 
0.541518926620483, 0.597665846347809, 0.320436120033264, 0.31931933760643, 
0.311251223087311, 0.311251223087311, 0.311251223087311, 0.317184776067734, 
0.338966459035873, 0.375375002622604, 0.367403119802475, 0.320436120033264, 
0.320436120033264, 0.488717943429947, 0.311251223087311, 0.311251223087311, 
0.427019357681274, 0.320436120033264, 0.316163510084152, 0.408030122518539, 
0.676600515842438, 0.65798020362854, 0.663134813308716, 0.405547618865967, 
0.666221380233765, 0.584704995155334, 0.519161760807037, 0.679032862186432, 
0.663134813308716, 0.493949443101883, 0.520084738731384, 0.519354522228241, 
0.673770666122437, 0.541518926620483, 0.592103779315948, 0.320436120033264, 
0.631685137748718, 0.663134813308716, 0.673770666122437, 0.573903799057007, 
0.38736030459404, 0.475033223628998, 0.663134813308716, 0.608104407787323, 
0.679032862186432, 0.657724738121033, 0.596750199794769, 0.634064376354218, 
0.32214692234993, 0.679032862186432, 0.609701991081238, 0.663134813308716, 
0.663134813308716, 0.663134813308716, 0.679032862186432, 0.666221380233765, 
0.526929080486298, 0.663134813308716, 0.663134813308716, 0.663134813308716, 
0.357654452323914, 0.539961099624634, 0.679032862186432, 0.553646564483643, 
0.611478388309479, 0.639116942882538, 0.663134813308716, 0.663134813308716, 
0.679032862186432, 0.632321059703827, 0.679032862186432, 0.519354522228241
), age = c(61L, 29L, 32L, 68L, 66L, 39L, 36L, 39L, 30L, 33L, 
75L, 44L, 63L, 66L, 67L, 31L, 52L, 45L, 38L, 33L, 63L, 46L, 69L, 
62L, 64L, 33L, 44L, 53L, 57L, 60L, 42L, 67L, 36L, 68L, 66L, 70L, 
42L, 39L, 43L, 64L, 59L, 34L, 73L, 65L, 79L, 19L, 51L, 65L, 70L, 
71L, 64L, 87L, 64L, 69L, 74L, 80L, 65L, 65L, 77L, 75L, 77L, 58L, 
54L, 57L, 81L, 53L, 85L, 73L, 62L, 57L, 52L, 82L, 71L, 78L, 74L, 
91L, 67L, 62L, 80L, 63L, 82L, 64L, 59L, 60L, 68L, 62L, 65L, 82L, 
76L, 68L, 71L, 62L, 74L, 63L), daysInHospital = c(15L, 17L, 14L, 
12L, 19L, 15L, 17L, 5L, 5L, 4L, 16L, 15L, 16L, 19L, 20L, 18L, 
15L, 17L, 14L, 12L, 12L, 21L, 13L, 11L, 13L, 16L, 15L, 18L, 18L, 
16L, 15L, 18L, 17L, 9L, 8L, 15L, 19L, 15L, 29L, 5L, 5L, 0L, 7L, 
2L, 2L, 3L, 4L, 29L, 14L, 7L, 4L, 9L, 18L, 10L, 9L, 13L, 14L, 
5L, 5L, 9L, 5L, 6L, 13L, 9L, 7L, 6L, 6L, 1L, 8L, 6L, 5L, 1L, 
2L, 4L, 9L, 5L, 4L, 1L, 28L, 9L, 6L, 6L, 17L, 11L, 2L, 29L, 29L, 
2L, 1L, 0L, 5L, 10L, 0L, 1L)), row.names = c(NA, -94L), class = c("data.table", 
"data.frame"))

编辑：

运行N 和K 的另一个随机样本似乎可以工作：

[1] "Current k:  10 Current n:  80"
[1] "Current MC run:  1 out of:  5"
[1] "Current MC run:  2 out of:  5"
[1] "Current MC run:  3 out of:  5"
[1] "Current MC run:  4 out of:  5"
[1] "Current MC run:  5 out of:  5"
[1] "Current k:  42 Current n:  91"
[1] "Current MC run:  1 out of:  5"
[1] "Current MC run:  2 out of:  5"
[1] "Current MC run:  3 out of:  5"
[1] "Current MC run:  4 out of:  5"
[1] "Current MC run:  5 out of:  5"
[1] "Current k:  17 Current n:  27"
[1] "Current MC run:  1 out of:  5"
[1] "Current MC run:  2 out of:  5"
[1] "Current MC run:  3 out of:  5"
[1] "Current MC run:  4 out of:  5"
[1] "Current MC run:  5 out of:  5"
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

我不确定为什么最后会得到 3 个 NULL。

【问题讨论】：

标签： r

【解决方案1】：

在第一个示例中出现错误的原因是因为您采样的 k 值大于 df (k = 96) 的行。因此，当您使用 k 的值对 df 的行进行采样时，它会出错。在您的第二个示例中不会出现此错误，因为k 和n 的所有值都小于94。

您有两种解决方案（不确定哪一种在您的应用中更好）。

在sample_n 函数中设置replace = TRUE。这样，您将进行替换抽样（但请注意，可能存在重复值）。
创建K 和N 的值时，在runif 函数中将上限设置为94 或更小。

请注意，当您使用slice(n) 时，这也可能是一个问题。所以可能第二种解决方案会更好。

您在第二个示例末尾获得 3 个 NULL 值的原因是因为 foreach 是一个返回最后执行命令的值的函数。在您的情况下，最后执行的命令位于 for 循环内，并且不会返回。所以foreach 返回 NULL。要隐藏这些返回值，您可以将foreach 的返回值放入一个虚拟变量中。

dummy.variable <- foreach(k = K, n = N) %do% {
  print(paste("Current k: ", k, "Current n: ", n))
  for(indMC in 1:MCtot){
    print(paste("Current MC run: ", indMC, "out of: ", MCtot ))
  }
}

只是为了说明foreach 的行为。最后一个命令返回 1 的示例。

foreach(k = K, n = N) %do% {
  print(paste("Current k: ", k, "Current n: ", n))
  for(indMC in 1:MCtot){
    print(paste("Current MC run: ", indMC, "out of: ", MCtot ))
  }
  1
}

[1] "Current k:  31 Current n:  42"
[1] "Current MC run:  1 out of:  5"
[1] "Current MC run:  2 out of:  5"
[1] "Current MC run:  3 out of:  5"
[1] "Current MC run:  4 out of:  5"
[1] "Current MC run:  5 out of:  5"
[1] "Current k:  61 Current n:  74"
[1] "Current MC run:  1 out of:  5"
[1] "Current MC run:  2 out of:  5"
[1] "Current MC run:  3 out of:  5"
[1] "Current MC run:  4 out of:  5"
[1] "Current MC run:  5 out of:  5"
[1] "Current k:  58 Current n:  81"
[1] "Current MC run:  1 out of:  5"
[1] "Current MC run:  2 out of:  5"
[1] "Current MC run:  3 out of:  5"
[1] "Current MC run:  4 out of:  5"
[1] "Current MC run:  5 out of:  5"
[[1]]
[1] 1

[[2]]
[1] 1

[[3]]
[1] 1

如果 1 值在内部 for 循环内，它将返回 NULL。

【讨论】：