Statistical sampling from a vector I have written a function/script to simulate data for sample size estimation. The following code samples from a vector of generated values with varying numbers of sample sizes and then concatenates means and standard deviations for a number of simulations.

I am using nested loops but would like to vectorize the code for efficiency and am wondering if anyone could suggest some optimizations.

library(MCMCglmm)
library(tidyverse)

Est &lt;- function(n, mean, sd, lower, upper, samp_min, samp_max, samp_int, nsim){

        Data &lt;- round(rtnorm(n, mean, sd, lower, upper), digits = 0) # Create a vector to sample from

        Samp_size &lt;- seq(samp_min, samp_max, samp_int) # Create vector of sample sizes

        # Set up enpty results data frames
        Results_samp &lt;- data.frame()
        Results &lt;- data.frame()

        for(i in 1:nsim){ ## Loop through number of simulations

                for (j in seq_along(Samp_size)) { # Loop through sample sizes
                     SUS_Score &lt;- sample(Data, Samp_size[j], replace = TRUE)
                     Nsubj &lt;- Samp_size[j]
                     Mean &lt;- mean(SUS_Score, na.rm = TRUE)
                     SD &lt;- sd(SUS_Score, na.rm = TRUE)
                     Results_samp &lt;- rbind(Results_samp, 
                                      data.frame(
                                              Nsubj,
                                              Mean,
                                              SD))
                }
             Results &lt;- rbind(Results, Results_samp) 
             Results_samp &lt;- data.frame()
        }

       Results %&gt;%
               arrange(Nsubj)
}

Sims &lt;- Est(n = 1000, mean = 55, sd = 37, lower = 0, upper = 100, samp_min = 5, samp_max = 35, samp_int = 5, nsim = 1000)


Any suggestions would be greatly appreciated!
Your code is quite good. Simulations is hard to vectorize. The largest slowdown here is the repeated calling of rbind in loop. It is faster to crate list of vectors and concatenate the results at the end. So  I edited the necessary parts:

Est3 &lt;- function(n, mean, sd, lower, upper, samp_min, samp_max, samp_int, nsim) {
  Data &lt;- round(rtnorm(n, mean, sd, lower, upper), digits = 0)
  Samp_size &lt;- seq(samp_min, samp_max, samp_int)
  Results &lt;- list() # crate emty list
  for (i in 1:nsim) {
    Results_samp &lt;- list() # crate emty list
    for (j in seq_along(Samp_size)) {
      Nsubj &lt;- Samp_size[j]
      SUS_Score &lt;- sample(Data, Nsubj, replace = TRUE)
      Mean &lt;- mean(SUS_Score, na.rm = TRUE)
      SD &lt;- sd(SUS_Score, na.rm = TRUE)
      Results_samp[[j]] &lt;- c(Nsubj, Mean, SD) # add values to list
    }
    Results[[i]] &lt;- Reduce(rbind, Results_samp) # convert list to matrix and add to main list
  }
  Results &lt;- Reduce(rbind, Results) # 'rbind' list of matrices
  Results &lt;- as.data.frame(Results, row.names = F) # convert to data.frame
  colnames(Results) &lt;- c('Nsubj', 'Mean', 'SD') # add names
  Results %&gt;% arrange(Nsubj)
}