How could I optimize this following program in R to boost performance? (Monte Carlo simulation involving compute-intensive permutation test)

In R, loops such as your for (i in 1:B) are slow. It's almost always best to avoid them. Here's a solution that avoids the loop entirely to perform a single permutation test.

start <- Sys.time()
permTest <- function(x) {
  # Create a tibble containing B replications of the sampled data
  tibble() %>% 
    expand(
      Perm=1:B,
      V=N$v
    ) %>% 
    # Randomly permute each permutation and hence randomly assign group membership 
    mutate(Random=runif(nrow(N) * B)) %>% 
    arrange(Perm, Random) %>% 
    mutate(g=rep(c("x", "y"), each=10, times=B)) %>% 
    # Calculate groupo means for each permutation 
    group_by(Perm, g) %>% 
    summarise(Mean=mean(V), .groups="drop") %>% 
    # Calculate p-value for each permutation and determine (one-sided) significance level
    pivot_wider(values_from=Mean, names_from=g) %>% 
    mutate(
      Delta=x-y,
      MoreExtreme=ifelse(Delta > D_obs, 1, 0)
    ) %>% 
    summarise(PValue=mean(MoreExtreme), .groups="drop") %>% 
    mutate(
      Power=ifelse(PValue < 0.05, 1, 0),
      Permutation=x
    )
}

This runs in 0.013s on my machine, suggesting that a thousand calls should take less than 15 seconds, approximately 30 times faster than the original. However,

# Simulate 1000 permutation tests
Result <- bind_rows(lapply(1:1000, permTest))

Took 3.38 minutes (~200s) suggesting that lapply and/or bind_rows are also slow. (Or that there was significant latency on my machine.) I expected some overhead dur to lapply and bind_rows, but I'm very surprised that the overhead is as large as it is. I don't have time to investigate further.

It is also possible to introduce the replication into the function directly:

permTest2 <- function(NTests) {
  # Create a tibble containing B replications of the sampled data
  tibble() %>% 
    expand(
      Perm=1:B,
      V=N$v,
      Test=1:NTests
    ) %>% 
    # Randomly permute each permutation and hence randomly assign group membership 
    mutate(Random=runif(nrow(N) * B)) %>% 
    arrange(Test, Perm, Random) %>% 
    mutate(g=rep(c("x", "y"), each=10, times=B)) %>% 
    # Calculate groupo means for each permutation 
    group_by(Test, Perm, g) %>% 
    summarise(Mean=mean(V), .groups="drop") %>% 
    # Calculate p-value for each permutation and determine (one-sided) significance level
    pivot_wider(values_from=Mean, names_from=g) %>% 
    mutate(
      Delta=x-y,
      MoreExtreme=ifelse(Delta > D_obs, 1, 0)
    ) %>% 
    group_by(Test) %>%
    summarise(PValue=mean(MoreExtreme), .groups="drop") %>% 
    mutate(
      Power=ifelse(PValue < 0.05, 1, 0)
    )
}

[Untested code]

But - on my machine at least - permTest2(1000) runs into memory problems. Therefore, the optimal soltion will likely be a combination of lapply and permTest2 with N less than 1000. In other words, use lapply to run chunks of the required permutations and then aggregate to get an overall result.

@denis' solution is also reasonable.

PS: your original post wasn't quite reproducible, but your edit whilst I was witing this answer may have solved that problem.

perm_test2 <- function(delta) {
  x <- rnorm(10, mean = delta, sd = 1)
  y <- rnorm(10, mean = 0, sd = 1)
  D_obs <- mean(x) - mean(y)
  v <- c(x, y)
  g <- rep(1:2, each = 10)
  B <- 5000
  y_mean <- x_mean <- 0
  D_perm <- rep(0, B)
  for (i in 1:B) {
    ii <- sample(g) == 1L
    D_perm[i] <- (sum(v[ii]) - sum(v[!ii]))/10 
  }
  P_perm <- sum(D_perm >= D_obs)/ B
  return(round(P_perm, digits = 5))
}

In comparison to more complicated approaches I propose the above one, where I have made some simple improvements to your existing code.

1. we do not need to use data.frame, that only takes unnecessary time to subset vector each time
2. instead of factor group vector we can use simple integer vector
3. the ordering of data is not needed
4. we can reduce mean calculation. In your code 2 calls of tapply(N$v, N$g, mean) was the slowest part.
5. mean(x) is slower than sum(x)/n, it does additional checks etc., so in this situation we can use the faster approach. As the inner loop will be executed 1000x5000 times(sims x B).

bench::mark(perm_test_org(0.5), perm_test2(0.5), iterations = 5, check = F,
            relative = T)[, 1:8]
# A tibble: 2 x 8
#   expression           min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
#   <bch:expr>         <dbl>  <dbl>     <dbl>     <dbl>    <dbl> <int> <dbl>
# 1 perm_test_org(0.5)  15.5   15.3       1        1.00     1        5    16
# 2 perm_test2(0.5)      1      1        13.5      1        1.69     5     2

Approximately 15x faster on my system. 1000 iterations took 33.89 seconds.

Update 2.

We can improve speed even more by:

1. replacing sample with sample.int
2. then we see that the g vector isn't needed at all to select random two groups
3. in loop we do not need to sum both parts of vector v as we can sum(v) before loop, so we can do one sum less inside the loop and calculate the result values at end.

perm_test3 <- function(delta) {
  x <- rnorm(10, mean = delta, sd = 1)
  y <- rnorm(10, mean = 0, sd = 1)
  D_obs <- mean(x) - mean(y)
  v <- c(x, y)
  B <- 5000
  s <- sum(v)
  D_perm2 <- rep(0, B)
  for (i in 1:B) {
    D_perm2[i] <- sum(v[sample.int(10) < 6])
  }
  D_perm <- D_perm2 - (s - D_perm2)
  P_perm <- sum(D_perm/10 >= D_obs) / B
  return(round(P_perm, digits = 5))
}

Runs in +/- 20 seconds for 1000 iterations. Now the slowest part is repeated call of sample.int. You can look into faster functions: https://www.r-bloggers.com/2019/04/fast-sampling-support-in-dqrng/

I have the following possibility, using data.table. The idea is to create one unique long data.table using CJ to then do the calculations in grouping operations.

power_p2  = function(numtrial,delta){
  x <- rnorm(10, mean = delta, sd = 1)
  y <- rnorm(10, mean = 0, sd = 1)
  D_obs <- mean(x) - mean(y)
  
  
  B <- 5000
  grouplist <- c(rep("x", 10), rep("y", 10))

  
  plouf <- CJ(it = 1:B,trial = 1:numtrial)
  Ntrial <- plouf[,.(v = c(x, y),g = sample(grouplist)),by = .(it,trial)]
  
  res <- Ntrial[,mean(v),by = .(g,it,trial)] %>%
    dcast(it + trial~g,value.var = "V1")
  pwr_p <- res[,.(p = sum(x-y>=D_obs)/B),by = trial][,sum(p <0.05)/numtrial]
  
  return(pwr_p)
}

The main gain of time is the use of CJ to generate the list of iteration from your first function, and trial from your second function.

plouf is thus a Bxnumtrial table. plouf[,.(v = c(x, y),g = sample(grouplist)),by = .(it,trial)] produce a data.table with the v = c(x, y) and the sampled group g = sample(grouplist) for each combination of it and trial.

The table Ntrial contains already the Bxnumtrial version of your sampling. The dcast is just to have a x and y column to easily calculate the difference of mean.

It is ~18 time faster, but I did not have the time to check properly that it does the same things than your function (I think it does, but I could have made an error)

microbenchmark(power_p(10,0.1),power_p2(10,0.1),times= 10 )

Unit: milliseconds
              expr        min         lq       mean     median         uq        max neval cld
  power_p(10, 0.1) 10873.7651 11524.0682 12240.7713 11748.8152 12389.2195 16029.2265    10   b
 power_p2(10, 0.1)   568.2784   600.1486   677.1702   632.8234   741.2503   992.4752    10  a