CODEWITHSUNDEEP
rfor-loopvectorization
KKW
KKW

Reputation: 385

how to write a more efficient code for this nested for-loop in R?

Trying to calculate a combination of p1 (event in experiment group) and p0 (event in control group) with or (odds ratio) of 1.5. nnt = Number needed to treat (100/(p1-p0))

library(tidyverse)

p1 <- seq(0,1, 0.0001)
p0 <- seq(0,1,0.0001)
or <- 1.5

df <- tibble(p1 = as.numeric(), p0 = as.numeric(), nnt = as.numeric())

for (i in p1) {
for (j in p0) {
  or_formula <- round((i/(1-i))/(j/(1-j)),3)
  
  if (or_formula == or & !is.na(or_formula)) {
    df <- df %>% add_row(p1 = i, p0 = j, nnt = round(1/(i-j), digits = 0))
  }
  
}
}

Upvotes: 3

Views: 77

Answers (3)

jpdugo17
jpdugo17

Reputation: 7106

A dplyr/purrr approach:

library(tidyverse)

p1 <- seq(0,1, 0.0001)
p0 <- seq(0,1,0.0001)
or <- 1.5


test <- map_df(p1, function(p1_element){
            map( p1_element,
                 ~list( 
                     p_0 = p0,
                     p_1 = p1_element,
                     or_formula = (. / (1-.)) / (p0 / (1 - p0)) ) ) }) %>%
        mutate(or_formula = round(or_formula, 3)) %>%
        filter(!is.na(or_formula) & or_formula == or) %>% 
        mutate(nnt = round(1 / (p_1 - p_0), digits =  0))

test %>% select(-or_formula)    
 

## 6.531 sec elapsed 
## original for loop : 131.454 sec elapsed

Upvotes: 1

ThomasIsCoding
ThomasIsCoding

Reputation: 101247

Here is another base R option with expand.grid + subset (but not as fast as @akurn's outer solution)

na.omit(
  subset(
    transform(
      expand.grid(p1 = p1, p0 = p0),
      nnt = round(1 / (p1 - p0), 0),
    ),
    round((p1 / (1 - p1)) / (p0 / (1 - p0)), 3) == or
  )
)

Upvotes: 1

akrun
akrun

Reputation: 887048

We could use outer

or_formula <- function(i, j) round((i/(1-i))/(j/(1-j)), 3)
m1 <- outer(p1, p0, FUN = or_formula)
dim(m1)
#[1] 10001 10001
i1 <- m1 == or & !is.na(m1)
i2 <- which(i1, arr.ind = TRUE)
p1new <- p1[i2[,1]]
p0new <- p0[i2[,2]]
df1 <- tibble(p1 = p1new, p0 = p0new, nnt = round(1/(p1new-p0new), digits = 0))

Benchmarks

-using outer

system.time({
  m1 <- outer(p1, p0, FUN = or_formula)
  i1 <- m1 == or & !is.na(m1)
  i2 <- which(i1, arr.ind = TRUE)
  p1new <- p1[i2[,1]]
  p0new <- p0[i2[,2]]
  df1 <- tibble(p1 = p1new, p0 = p0new, nnt = round(1/(p1new-p0new), digits = 0))
 
 
 })
#   user  system elapsed 
#  5.038   1.288   6.319

-using OP's for loop

system.time({
  df <- tibble(p1 = as.numeric(), p0 = as.numeric(), nnt = as.numeric())
  
  for (i in p1) {
   for (j in p0) {
   or_formula <- round((i/(1-i))/(j/(1-j)),3)
   
   if (or_formula == or & !is.na(or_formula)) {
     df <- df %>% add_row(p1 = i, p0 = j, nnt = round(1/(i-j), digits = 0))
   }
   
 }
 }
 
 
 })
#   user  system elapsed 
#122.391   0.748 123.128

-testing the equality

identical(df, df1)
#[1] TRUE

Upvotes: 4

Related Questions