Using parallelization

Simple parallelization without custom functions

If we use standard model fitting function (e.g., “lm”) that are included in the base package, parallelization is comparatively simple. We only need to load the package furrr and specify a “plan for how to resolve a future” (for more information see ?future::plan). In this case, I am choosing multisession (resolve the computations separate R sessions running in the background on the same machine) and specify workers = 4 so that it runs on 4 cores in parallel.

# Setup of specifications (number of specs = 1152)
specs <- setup(data = dat,            
               y = c("y1", "y2", "y3"),                    
               x = c("x1", "x2", "x3", "x4"),               
               model = c("lm"), 
               controls = c("c1", "c2", "c3", "c4"),
               subsets = list(group = unique(dat$group)))   

# Default: Sequential ---
results_simple <- specr(specs)

# Parallel: Multisession (only works when `furrr` is loaded!)
plan(strategy = multisession, workers = 4)
results_parall <- specr(specs)

# Comparison
cat("Sequential: ", results_simple$time, "\n",
    "Parallel:   ", results_parall$time)

## Sequential:  57.227 sec elapsed 
## Parallel:    36.781 sec elapsed

As we can see, the default (sequential) computation took around a minute (57.227 sec elapsed) and the parallel computation about half a minute (36.781 sec elapsed).

We have to acknowledge that even with this comparatively large data set and more than 1,000 specifications, the reduction in time not too much. Thus, parallelization only makes sense if you have a truly large data set and thousands of specifications or the type of model you are estimating takes a really long time (e.g., a complex structural equation model, a negative binomial model, etc.). The true power of parallelization furthermore only come into play if you are using a considerable number of cores.

Parallelization with custom functions from different packages

If we use a custom function, we need to make sure that this function is passed to the different workers. This can be done by specifying so called furrr_options. We need to pass objects from the global environment (and potentially also packages) as shown below. Please note that we do not have to specify the “future plan” again, because have specified it already earlier in this session (see above). If we are unsure what plan is currently specified, we can simply run plan() and get some information about the current setup.

As computation can take a long time, it would be nice if we would see some kind of progress indication. This can easily be done by simply adding the argument .progress = TRUE to specr(). This passes this argument to the future_pmap() function within specr() and prints a rudimentary progress bar during the fitting process.

# Custom function
log_model <- function(formula, data) {
  glm(formula, data, family = binomial())
}

# Setup specs
specs <- setup(data = dat,            
               y = c("y4"),                    
               x = c("x1", "x2", "x3"),               
               model = c("log_model"), 
               controls = c("c1", "c2"),
               subsets = list(group = unique(dat$group)))   

# Create furrr_options to be passed to specr() (only works if `furrr` is loaded)
opts <- furrr_options(
  globals = list(log_model = log_model)
)

# What "plan" is currently specified?
plan()

## multisession:
## - args: function (..., workers = 4, envir = parent.frame())
## - tweaked: TRUE
## - call: plan(strategy = multisession, workers = 4)

# Run results
results_parall_2 <- specr(specs, 
                          .options = opts,   # Pass ops to specr
                          .progress = TRUE)  # To add progress bar (not shown in output)

# Summarize results
summary(results_parall_2)

## Results of the specification curve analysis
## -------------------
## Technical details:
## 
##   Class:                          specr.object -- version: 1.0.0 
##   Cores used:                     4 
##   Duration of fitting process:    37.736 sec elapsed 
##   Number of specifications:       72 
## 
## Descriptive summary of the specification curve:
## 
##  median  mad   min  max q25  q75
##    0.01 0.01 -0.02 0.03   0 0.01
## 
## Descriptive summary of sample sizes: 
## 
##  median   min   max
##   20038 19844 1e+05
## 
## Head of the specification results (first 6 rows): 
## 
## # A tibble: 6 × 21
##   x     y     model  controls subsets group formula estimate std.error statistic p.value conf.low 
##   <chr> <chr> <chr>  <chr>    <chr>   <fct> <glue>     <dbl>     <dbl>     <dbl>   <dbl>    <dbl> 
## 1 x1    y4    log_m… no cova… e       e     y4 ~ x…     0.01      0.01      0.38    0.7     -0.02 
## 2 x1    y4    log_m… no cova… a       a     y4 ~ x…     0.01      0.01      0.62    0.54    -0.02 
## 3 x1    y4    log_m… no cova… d       d     y4 ~ x…    -0.02      0.01     -1.17    0.24    -0.04 
## 4 x1    y4    log_m… no cova… c       c     y4 ~ x…     0         0.01     -0.16    0.87    -0.03 
## 5 x1    y4    log_m… no cova… b       b     y4 ~ x…     0.03      0.01      1.84    0.07     0    
## 6 x1    y4    log_m… no cova… all     <NA>  y4 ~ x…     0         0.01      0.66    0.51    -0.01 
## # … with 9 more variables: conf.high <dbl>, fit_null.devian… <dbl>, fit_df.null <dbl>, 
## #   fit_logLik <dbl>, fit_AIC <dbl>, fit_BIC <dbl>, fit_deviance <dbl>, fit_df.residual <dbl>, 
## #   fit_nobs <dbl>

Note: In the technical details of the summary, we also always see how many cores were used and how long the fitting process has taken.

At the end of our analysis, it makes sense to explicitly close multisession workers by switching the plan back to sequential.

plan(sequential)