How to iterate through parameters in for loop

Viewed 1351

I have a model written as a for loop that incorporates a number of parameters that I specify:

## functions needed to run the model
learn <- function(prior, sensi, speci, e){
  out <- ifelse(e == 1, (sensi*prior) / ((sensi*prior) + (1-speci)*(1-prior)),
                ((1-sensi)*prior) / (((1-sensi)*prior) + (speci*(1-prior))))
  out
}

feed <- function(vec){
  prior <- 0.5
  for (i in vec){
    res <- learn(prior, sensi, speci, i)
    prior <- res
  }
  return(prior)
}

## specify parameters
iterations <- 100
N <- 10
BR <- 0.66
sensi <- 0.75
speci <- 0.45

## initialize results object
res <- NULL

## loop for number of iterations
for (j in 1:iterations){
  
  X <- as.numeric(rbinom(1, 1, BR))
  
  if (X == 1){ # if X is 1...
    agents <- c(1:N) 
    evidence <- vector("list", length(agents)) 
    for (i in agents) {
      n <- sample(10, 1, replace = TRUE) 
      evidence[[i]] <- rbinom(n, 1, sensi) 
    }
  } else { # if X is 0... 
    agents <- c(1:N)
    evidence <- vector("list", length(agents)) 
    for (i in agents) {
      n <- sample(10, 1, replace = TRUE) 
      evidence[[i]] <- rbinom(n, 1, sensi) 
      evidence[[i]] <- ifelse(evidence[[i]]==1, 0, 1) # flip evidence 
    }
  }
  
  # feed vectors of evidence through learn function
  t0 <- sapply(evidence, feed)
  
  # save dataframe 
  df <- data.frame("i" = j, 
                   "ID" = c(1:N), 
                   "E" = t0, 
                   "X" = X,
                   "N" = N, 
                   "BR" = BR,
                   "sensi" = sensi,
                   "speci" = speci)

  res <- rbind(res, df)
  
}

This works fine for a single parameterisation, but I now want to automate the process of specifying different parameter values and re-running the model. So instead of defining each parameter as a single value, I define them as a vector of values and store all the possible parameterisations in a dataframe (paramspace) with each row holding the values for a single parameterisation that I want to run:

## set up for multiple parameterizations 
iterations <- 100
N_vec <- c(10, 50)
BR_vec <- c(0.25, 0.50, 0.75) 
sensi_vec <- c(0.45, 0.75)
speci_vec <- c(0.45, 0.75)

paramspace <- expand.grid(iterations = iterations, N = N_vec, BR = BR_vec, sensi = sensi_vec, speci = speci_vec)

> paramspace
   iterations  N   BR sensi speci
1         100 10 0.25  0.45  0.45
2         100 50 0.25  0.45  0.45
3         100 10 0.50  0.45  0.45
4         100 50 0.50  0.45  0.45
5         100 10 0.75  0.45  0.45
6         100 50 0.75  0.45  0.45
7         100 10 0.25  0.75  0.45
8         100 50 0.25  0.75  0.45
9         100 10 0.50  0.75  0.45
10        100 50 0.50  0.75  0.45
11        100 10 0.75  0.75  0.45
12        100 50 0.75  0.75  0.45
13        100 10 0.25  0.45  0.75
14        100 50 0.25  0.45  0.75
15        100 10 0.50  0.45  0.75
16        100 50 0.50  0.45  0.75
17        100 10 0.75  0.45  0.75
18        100 50 0.75  0.45  0.75
19        100 10 0.25  0.75  0.75
20        100 50 0.25  0.75  0.75
21        100 10 0.50  0.75  0.75
22        100 50 0.50  0.75  0.75
23        100 10 0.75  0.75  0.75
24        100 50 0.75  0.75  0.75

How can I pass each row of parameter values to my model and automatically run through all the parameterisations stated in paramspace?

3 Answers

As suggested in comments, you can create a function and then use apply to loop over the parameters combinations :


## functions needed to run the model
learn <- function(prior, sensi, speci, e){
  out <- ifelse(e == 1, (sensi*prior) / ((sensi*prior) + (1-speci)*(1-prior)),
                ((1-sensi)*prior) / (((1-sensi)*prior) + (speci*(1-prior))))
  out
}

feed <- function(vec,sensi,speci){
  prior <- 0.5
  for (i in vec){
    res <- learn(prior, sensi, speci, i)
    prior <- res
  }
  return(prior)
}

runModel <- function(iterations = 100,
                     N = 10,
                     BR = 0.66,
                     sensi = 0.75,
                     speci = 0.45 ) {
  ## initialize results object
  res <- NULL
  
  ## loop for number of iterations
  for (j in 1:iterations){
    
    X <- as.numeric(rbinom(1, 1, BR))
    
    if (X == 1){ # if X is 1...
      agents <- c(1:N) 
      evidence <- vector("list", length(agents)) 
      for (i in agents) {
        n <- sample(10, 1, replace = TRUE) 
        evidence[[i]] <- rbinom(n, 1, sensi) 
      }
    } else { # if X is 0... 
      agents <- c(1:N)
      evidence <- vector("list", length(agents)) 
      for (i in agents) {
        n <- sample(10, 1, replace = TRUE) 
        evidence[[i]] <- rbinom(n, 1, sensi) 
        evidence[[i]] <- ifelse(evidence[[i]]==1, 0, 1) # flip evidence 
      }
    }
    
    # feed vectors of evidence through learn function
    #t0 <- sapply(evidence, feed)
    t0 <- sapply(evidence,function(e){feed(e,sensi,speci)})
    
    # save dataframe 
    df <- list("i" = iterations, 
               "ID" = c(1:N), 
               "E" = t0, 
               "X" = X,
               "N" = N, 
               "BR" = BR,
               "sensi" = sensi,
               "speci" = speci)
    
    res <- rbind(res, df)
    
  }
  res
}

# Define parameter space
iterations <- 100
N_vec <- c(10, 50)
BR_vec <- c(0.25, 0.50, 0.75) 
sensi_vec <- c(0.45, 0.75)
speci_vec <- c(0.45, 0.75)

paramspace <- expand.grid(iterations = iterations, N = N_vec, BR = BR_vec, sensi = sensi_vec, speci = speci_vec)

# Loop over parameter space :
res <- apply(paramspace,1,function(paramset) {
  iterations = paramset[1]
  N = paramset[2]
  BR = paramset[3]
  sensi = paramset[4]
  speci = paramset[5]
  runModel(iterations = iterations, N = N, BR = BR , sensi = sensi, speci = speci )
})

You can also use the foreach package, that used with an appropriate backend offers parallelization capabilities, in case your task becomes more intensive. Here a simple example to understand how it works.

foreach(a=1:3, b=4:6) %do% (a + b)

Then I tried to embed your code into foreach


require(foreach)

## functions needed to run the model
learn <- function(prior, sensi, speci, e){
  out <- ifelse(e == 1, (sensi*prior) / ((sensi*prior) + (1-speci)*(1-prior)),
                ((1-sensi)*prior) / (((1-sensi)*prior) + (speci*(1-prior))))
  out
}

feed <- function(vec){
  prior <- 0.5
  for (i in vec){
    res <- learn(prior, sensi, speci, i)
    prior <- res
  }
  return(prior)
}


## set up for multiple parameterizations 
iterations <- 100
N_vec <- c(10, 50)
BR_vec <- c(0.25, 0.50, 0.75) 
sensi_vec <- c(0.45, 0.75)
speci_vec <- c(0.45, 0.75)

paramspace <- expand.grid(iterations = iterations, N = N_vec, BR = BR_vec, sensi = sensi_vec, speci = speci_vec)

res <- foreach(iterations = paramspace$iterations, 
               N = paramspace$N, 
               BR = paramspace$BR, 
               sensi = paramspace$sensi, 
               speci = paramspace$speci) %do% {
                 
                 ## initialize results object
                 res <- NULL
                 
                 ## loop for number of iterations
                 for (j in 1:iterations){
                   
                   X <- as.numeric(rbinom(1, 1, BR))
                   
                   if (X == 1){ # if X is 1...
                     agents <- c(1:N) 
                     evidence <- vector("list", length(agents)) 
                     for (i in agents) {
                       n <- sample(10, 1, replace = TRUE) 
                       evidence[[i]] <- rbinom(n, 1, sensi) 
                     }
                   } else { # if X is 0... 
                     agents <- c(1:N)
                     evidence <- vector("list", length(agents)) 
                     for (i in agents) {
                       n <- sample(10, 1, replace = TRUE) 
                       evidence[[i]] <- rbinom(n, 1, sensi) 
                       evidence[[i]] <- ifelse(evidence[[i]]==1, 0, 1) # flip evidence 
                     }
                   }
                   
                   # feed vectors of evidence through learn function
                   t0 <- sapply(evidence, feed)
                   
                   # save dataframe 
                   df <- data.frame("i" = j, 
                                    "ID" = c(1:N), 
                                    "E" = t0, 
                                    "X" = X,
                                    "N" = N, 
                                    "BR" = BR,
                                    "sensi" = sensi,
                                    "speci" = speci)
                   
                   res <- rbind(res, df)
                   
                 }
                 
                 res
                 
               }

Another approach is to make a function and to use Map(...). The advantage of Map is that your paramspace will not be coerced into a matrix which will make everything the same type (i.e., numeric, character, etc.).

There were also some other changes I made in order to allow R to do the acccounting for us. Primarily:

  1. X is now a logical so we can simplify our if statements. Additionally, the allocation is made all at once instead of looping.
  2. We change the feed() function to also generate the evidence. This allows us to...
  3. Use replicate to repeat the loops.
learn2 <- function(prior, sensi, speci, e){
    out <- ifelse(e, (sensi*prior) / ((sensi*prior) + (1-speci)*(1-prior)),
                  ((1-sensi)*prior) / (((1-sensi)*prior) + (speci*(1-prior))))
    out
}


feed2  = function(x, N, samp_n = 10L, sensi, speci) {
    evidence = rbinom(sample(samp_n, 1L, replace = TRUE), 
                      1,
                      if (x) sensi else 1 - sensi)
    
    prior = 0.5
    for (i in evidence) {
        res = learn2(prior, sensi, speci, i)
        prior = res
    }
    return(prior)
}

runModel2 <- function(iterations = 2,
                     N = 10,
                     BR = 0.66,
                     sensi = 0.75,
                     speci = 0.45 ) {
    
    X = sample(c(TRUE, FALSE), N, BR)
    
    ## this is done now so that the columns will be ordered nicer
    ans = list(ID = 1:N,
                N = N,
               BR = BR,
               sensi = sensi,
               speci = speci,
               X = X)
    
    t0s = replicate(iterations,
                    vapply(X, feed2, FUN.VALUE = 0, N, 10L, sensi, speci, USE.NAMES = FALSE), 
                    simplify = FALSE)
    
    names(t0s) = paste0("E_", 1:iterations)
    
    return(as.data.frame(c(ans, t0s)))
}

runModel2()
#>    ID  N   BR sensi speci     X        E_1         E_2
#> 1   1 10 0.66  0.75  0.45  TRUE 0.82967106 0.657648599
#> 2   2 10 0.66  0.75  0.45 FALSE 0.43103448 0.006827641
#> 3   3 10 0.66  0.75  0.45  TRUE 0.43103448 0.775671866
#> 4   4 10 0.66  0.75  0.45  TRUE 0.71716957 0.431034483
#> 5   5 10 0.66  0.75  0.45 FALSE 0.24176079 0.016593958
#> 6   6 10 0.66  0.75  0.45 FALSE 0.30303324 0.008992838
#> 7   7 10 0.66  0.75  0.45  TRUE 0.82967106 0.865405260
#> 8   8 10 0.66  0.75  0.45 FALSE 0.43103448 0.439027817
#> 9   9 10 0.66  0.75  0.45 FALSE 0.57692308 0.050262167
#> 10 10 10 0.66  0.75  0.45 FALSE 0.02178833 0.296208531

This output is a little wider than your original approach. We can always reshape the E_# columns but this may end up being better for your actual use case.

Finally, here is Map() in action:

iterations <- 100
N_vec <- c(10, 50)
BR_vec <- c(0.25, 0.50, 0.75) 
sensi_vec <- c(0.45, 0.75)
speci_vec <- c(0.45, 0.75)

paramspace <- expand.grid(iterations = iterations, N = N_vec, BR = BR_vec, sensi = sensi_vec, speci = speci_vec)

res = Map(runModel2, paramspace$iterations, paramspace$N, paramspace$BR, paramspace$sensi, paramspace$speci)

res[[24L]][1:10, 1:8] ## only first 10 rows for demonstration
##   ID  N   BR sensi speci     X         E_1         E_2
##1   1 50 0.75  0.75  0.75  TRUE 0.500000000 0.500000000
##2   2 50 0.75  0.75  0.75 FALSE 0.001369863 0.035714286
##3   3 50 0.75  0.75  0.75 FALSE 0.250000000 0.900000000
##4   4 50 0.75  0.75  0.75  TRUE 0.750000000 0.250000000
##5   5 50 0.75  0.75  0.75  TRUE 0.987804878 0.500000000
##6   6 50 0.75  0.75  0.75  TRUE 0.964285714 0.250000000
##7   7 50 0.75  0.75  0.75  TRUE 0.750000000 0.750000000
##8   8 50 0.75  0.75  0.75 FALSE 0.012195122 0.035714286
##9   9 50 0.75  0.75  0.75  TRUE 0.750000000 0.500000000
##10 10 50 0.75  0.75  0.75 FALSE 0.250000000 0.001369863
Related