rperformanceparallel-processingregression# Does it make sense to use parallel with a bayesian regression?

I am runing a bayesian regression in R with 1000 chains and it is taking to long. In my attempt to improve the performance I tried parallel. I read that this technique is useful only when there are tasks that can be performed independently from each other. As I see the Markov Chains as iterations, I belived that this was the case. However, I am getting a worse speed compared to not using parallel.

If the problem is that parallel is not useful in this case. Is there another way to make it faster?

These are the libraries and replicable datasets:

```
library(rstan)
library(rstanarm)
library(tidyverse)
library(stringr)
library(openxlsx)
library(tictoc)
library(parallel)
tasas_ori <- data.frame(
gedad = c(rep(c("A", "B", "C", "D", "E", "F"), 6)),
employment = c(rep(c(20:25), 6),
rep(c(20:25) + 1, 6),
rep(c(20:25) + 2, 6),
rep(c(20:25) + 3, 6),
rep(c(20:25) + 4, 6),
rep(c(20:25) + 5, 6)),
anio = c(rep(2006, 6), rep(2012, 6), rep(2014, 6),
rep(2018, 6), rep("2018 - aj.", 6), rep(2023, 6))
)
tasas <- tasas_ori %>%
filter(!(anio %in% c("2018-aj.", "2018"))) %>%
mutate(anio = as.numeric(anio) - 2006)
tasas2 <- tasas_ori %>%
filter(!(anio %in% c("2018-aj."))) %>%
mutate(anio = as.numeric(anio) - 2006)
```

This is the code with parallel:

```
set.seed(1984)
tic()
num_cores <- parallel::detectCores()
cl <- makeCluster(num_cores - 2)
registerDoParallel(cl)
for (i in 1:3) {
# Parallel loop to fit stan_glm models
bayes_reg <- foreach(i = 1:1, .combine = "c") %dopar% {
library(rstan) # Load rstan within the parallel worker
library(rstanarm) # Load rstanarm within the parallel worker
library(dplyr)
stan_glm(as.formula(paste("employment", "~ anio")),
data = tasas %>% filter(gedad == unique(tasas$gedad)[1]),
chain = 50,
# prior = normal(normal_reg$coefficients[2], summary(normal_reg)$coefficients[4]),
refresh = 0
)
}
# Stop the parallel processing
summary(bayes_reg, digits = 5) %>% print()
coefficients_by[i] <- as.data.frame(bayes_reg$coefficients)
predictions_by[i] <- as.data.frame(predict(bayes_reg, newdata = tasas2 %>%
filter(gedad == unique(tasas$gedad)[i])))
}
stopCluster(cl)
toc()
```

And this is the code without parallel:

```
# Regresión por grupos etarios
coefficients_by <- vector("list")
predictions_by <- vector("list")
bayes <- vector("list")
set.seed(1984)
tic()
for (i in 1:3) {
cat(paste("\n", "**** Regresión bayesiana", "-", "Grupo de edad: ", unique(dci_tasas$gedad)[i], "****", "\n"))
bayes_reg <- stan_glm(as.formula(paste("employment", "~ anio")),
data = dci_tasas %>%
filter(gedad == unique(dci_tasas$gedad)[i]),
chain = 50,
# prior_intercept = normal(normal_reg$coefficients[1], summary(normal_reg)$coefficients[3]),
# prior = normal(normal_reg$coefficients[2], summary(normal_reg)$coefficients[4]),
refresh = 0
)
summary(bayes_reg, digits = 5) %>% print()
coefficients_by[i] <- as.data.frame(bayes_reg$coefficients)
predictions_by[i] <- as.data.frame(predict(bayes_reg, newdata = dci_tasas2 %>%
filter(gedad == unique(dci_tasas$gedad)[i])))
}
toc()
stopCluster(cl)
toc()
```

Solution

I think the problem is that you're not actually using `foreach()`

as it's intended. The `foreach()`

function should execute the loop. Looping over 1 value really means that it's not running in parallel. In fact, you're suffering all the startup time of making a cluster without generating any of the benefits. If you use `foreach()`

to do the looping, you'll see you get not quite a 100% increase in efficiency.

```
library(rstan)
library(rstanarm)
library(dplyr)
library(tictoc)
tasas_ori <- data.frame(
gedad = c(rep(c("A", "B", "C", "D", "E", "F"), 6)),
employment = c(rep(c(20:25), 6),
rep(c(20:25) + 1, 6),
rep(c(20:25) + 2, 6),
rep(c(20:25) + 3, 6),
rep(c(20:25) + 4, 6),
rep(c(20:25) + 5, 6)),
anio = c(rep(2006, 6), rep(2012, 6), rep(2014, 6),
rep(2018, 6), rep(2019, 6), rep(2023, 6))
)
tasas <- tasas_ori %>%
filter(!(anio %in% c(2018, 2019))) %>%
mutate(anio = as.numeric(anio) - 2006)
```

In the parallel code below, I save each result as a list and then the output from `foreach(i=1:3, ...) %dopar% ...`

is a list with three elements - one of reach of the models run. I setup a cluster with three cores (since the loop is over three values, more cores wouldn't matter). The result takes 4.6 seconds to computer.

```
set.seed(1984)
library(foreach)
library(doParallel)
#> Loading required package: iterators
#> Loading required package: parallel
pcl<-makeCluster(3)
registerDoParallel(pcl)
tic()
bayes_reg <- foreach(i = 1:3, .combine = "c") %dopar% {
library(rstanarm)
library(dplyr)
tasas_ori <- data.frame(
gedad = c(rep(c("A", "B", "C", "D", "E", "F"), 6)),
employment = c(rep(c(20:25), 6),
rep(c(20:25) + 1, 6),
rep(c(20:25) + 2, 6),
rep(c(20:25) + 3, 6),
rep(c(20:25) + 4, 6),
rep(c(20:25) + 5, 6)),
anio = c(rep(2006, 6), rep(2012, 6), rep(2014, 6),
rep(2018, 6), rep(2019, 6), rep(2023, 6))
)
tasas <- tasas_ori %>%
filter(!(anio %in% c(2018, 2019))) %>%
mutate(anio = as.numeric(anio) - 2006)
res <- stan_glm(as.formula(paste("employment", "~ anio")),
data = tasas %>% filter(gedad == unique(tasas$gedad)[1]),
chain = 50,
# prior = normal(normal_reg$coefficients[2], summary(normal_reg)$coefficients[4]),
refresh = 0
)
list(res)
}
stopCluster(pcl)
toc()
#> 4.62 sec elapsed
```

In the code below, I use a `for()`

loop without parallelization, just like you did. You can see that the result takes 7.4 seconds to compute. So, the for loop is 7.4/4.6 = 1.6 or 60% slower. This is about what you might expect. The setup and communication overhead means that you're not going to get 3x the performance for 3x the cores.

```
tic()
bayes_reg2 <- vector(mode="list", length=3)
for(i in 1:3){
bayes_reg2[[i]] <- stan_glm(as.formula(paste("employment", "~ anio")),
data = tasas %>% filter(gedad == unique(tasas$gedad)[1]),
chain = 50,
# prior = normal(normal_reg$coefficients[2], summary(normal_reg)$coefficients[4]),
refresh = 0)
}
toc()
#> 7.386 sec elapsed
```

^{Created on 2023-09-03 with reprex v2.0.2}

- Installing R on Linux: configure: error: libcurl >= 7.28.0 library and headers are required with support for https
- How to do ensembles with time series using AICc?
- planes3d expands and draws the area based on the sphere's radius
- How to extract tag code itself using R, rvest
- How to Display or Print Contents of Environment in R
- How to use Windows user credentials for proxy authentication in R/RStudio
- R reticulate specifying python executable to use
- Replace multiple Instances of a variable name in an R function and save the modified function
- Standardizing address formatting in R
- How to fix "failed to load cairo DLL" in R?
- Using grepl to filter columns names in specific range of columns
- changing the legends in ggplot2 to have groups of similar labels
- How to keep only unique rows but ignore a column?
- convert string date to R Date FAST for all dates
- Add subgroup text to plotly pie chart
- R Shiny : adjust height of DT datatable when fillContainer=TRUE,
- Why do R external pointers' "unusual copying semantics" mean they should not be used stand-alone?
- How to extract somo character after a string with a number of word which can change in R
- What does `se` stand for in geom_smooth(..., se = FALSE)?
- How to find number of rows greater than any values in R
- Align text and reduce space between text and parentheses in plotly hover info box
- Remove outer box of geom_bar plot with broken y-axis
- How to use lag/lead in mutate with an initial value?
- Is it possible to have a Shiny ConditionalPanel whose condition is a global variable?
- counting elements in one list in another list
- How to vectorize nested loops in R?
- Replace NA values with an incrementing sequence starting from the previous non-NA value
- How can I calculate the number of uniques in a row within a species matrix?
- How to perform operations on pairs of rows, based on a "distinguishing" column's values
- Mutate variable based on previous observations