# Running compiled Stan models in Shiny

**R-bloggers | A Random Walk**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# Introduction

The aim of this post is to provide a short step-by-step guide on writing interactive R Shiny-applications that include models written in Stan using `rstan`

and `rstantools`

. The remainder of this post assumes a small amount of working knowledge on writing models in Stan and usage of the package `rstan`

to interface Stan from R.

## Demo set-up

For demonstration purposes, let’s start by writing a minimal Stan model file `lm.stan`

:

data { int<lower=0> N; vector[N] x; vector[N] y; } parameters { real alpha; real beta; real<lower=0> sigma; } model { y ~ normal(alpha + beta * x, sigma); }

This Stan file encodes a simple linear regression model of the form:

\[ y_i \ \overset{\text{iid}}{\sim} \ N(\alpha + \beta \cdot x_i, \sigma^2), \quad i = 1,\ldots,N \]

Next, we create a small shiny-app contained in a single file `app.R`

(in the same directory as `lm.stan`

) that draws posterior samples from the Stan model in `lm.stan`

with calls to `rstan::stan_model()`

and `rstan::sampling()`

:

library(shiny) library(rstan) ui <- fluidPage( sidebarLayout( sidebarPanel( numericInput("N", label = "N", value = 10) ), mainPanel( plotOutput("posteriors") ) ) ) server <- function(input, output, session) { ## compile stan model model <- stan_model(file = "lm.stan") ## draw samples draws <- reactive({ N <- input$N sampling( object = model, data = list(N = N, x = seq_len(N), y = rnorm(N, seq_len(N), 0.1)), chains = 2, iter = 1000 ) }) ## plot histograms output$posteriors <- renderPlot({ req(draws()) op <- par(mfrow = c(1, 2), cex = 1.25) hist(extract(draws(), "alpha")[[1]], main = bquote("Posterior samples"~alpha), xlab = expression(alpha)) hist(extract(draws(), "beta")[[1]], main = bquote("Posterior samples"~beta), xlab = expression(beta)) par(op) }) } shinyApp(ui = ui, server = server)

The contents of this shiny-app can be summarized in a simple reactive graph:

New noisy responses \(y_1, \ldots, y_N\) are generated according to \(y_i \overset{\text{iid}}{\sim} N(\alpha + \beta \cdot x_i, \sigma^2)\) with \(x_i = i\) for each \(i = 1,\ldots, N\) and fixed parameters \(\alpha = 0\), \(\beta = 1\) and \(\sigma = 0.1\).

## Slow model compilation

A fixed number of 2000 posterior samples for \(\alpha\) and \(\beta\) (and \(\sigma\)) is drawn across two individual chains (i.e. 1000 draws per chain), which should not take more than a few seconds to complete on an ordinary laptop computer, especially if \(N\) is small. However, when launching the shiny-app, it becomes evident that it takes significantly longer to complete drawing any initial posterior samples.

The reason for this is (obviously) that the Stan model has to be **recompiled** from the `lm.stan`

file whenever we launch the shiny-app in a new R-session due to the call to `rstan::stan_model()`

. Depending on the compiler settings, it takes up to ~1 minute to compile this single Stan model on my laptop computer, which more or less defeats the purpose of creating a shiny-app for interactive use.

Luckily, it is quite simple to avoid this unnecessary computational effort: we just have to compile our Stan models beforehand so that we can sample directly from the already compiled Stan models and skip the compilation step in `rstan::stan_model()`

. Before describing a general R-package approach using `rstantools`

, we start with a simpler approach –suitable for a large set of standard regression models– which is to take advantage of the pre-compiled Stan models in `rstanarm`

.

# Pre-compiled models with `rstanarm`

If the model we wish to sample from is already made available in R via the `rstanarm`

-package, arguably the most straightforward approach to avoid unnecessary Stan model compilation is to use `rstanarm`

’s R wrapper functions to directly sample from a pre-compiled Stan model. Note that if we are fitting a relative standard regression model, there is a decent chance a pre-compiled model version is available in `rstanarm`

. Besides the fact that the Stan models in `rstanarm`

are pre-compiled, the implementations of the Stan programs are likely more robust and computationally stable than any quick Stan program we would implement ourselves.

To sample from a simple linear model as defined in `lm.stan`

with `rstanarm`

, it suffices to remove the call to `rstan::stan_model()`

in `app.R`

and replace `rstan::sampling()`

by a call to `rstanarm::stan_lm()`

or `rstanarm::stan_glm()`

^{1}:

library(shiny) library(rstanarm) ui <- fluidPage( sidebarLayout( sidebarPanel( numericInput("N", label = "N", value = 10) ), mainPanel( plotOutput("posteriors") ) ) ) server <- function(input, output, session) { ## draw samples directly draws <- reactive({ N <- input$N samples <- stan_glm( formula = y ~ x, data = data.frame(x = seq_len(N), y = rnorm(N, seq_len(N), 0.1)), chains = 2, iter = 1000 ) as.matrix(samples)[, c(1, 2)] }) ## plot histograms output$posteriors <- renderPlot({ req(draws()) op <- par(mfrow = c(1, 2), cex = 1.25) hist(draws()[, 1], main = bquote("Posterior samples"~alpha), xlab = expression(alpha)) hist(draws()[, 2], main = bquote("Posterior samples"~beta), xlab = expression(beta)) par(op) }) } shinyApp(ui = ui, server = server)

The modified shiny-app no longer exhibits the same lack of responsiveness due to (unnecessary) Stan model recompilation, and it is seen from the animation that new posterior samples are generated almost instantly.

**Remark**: note that the model used by `stan_glm()`

is not exactly equivalent to the model in `lm.stan`

, since `stan_glm()`

assigns weakly informative priors to the model parameters by default^{2}, whereas (non-informative) uniform priors are used in the original `lm.stan`

file.

That’s great, but what if the model we wish to fit is not available through the `rstanarm`

-package? In this case, we can mimic the same general approach that `rstanarm`

follows: compile the Stan models **once** on package installation, and directly sample from the pre-compiled models thereafter. As it turns out, we can make use of the excellent `rstantools`

-package for exactly this purpose. The `rstantools`

-package essentially eliminates the effort of setting up a correct R-package structure and only requires us to include the Stan programs that should be compiled with the R-package.

# Creating a package with `rstantools`

To set up a new R-package that should interface Stan, we call `rstantools::rstan_create_package()`

which is roughly similar in use to `package.skeleton()`

, (or `usethis::create_package()`

or `Rcpp::Rcpp.package.skeleton()`

for that matter). The already existing `lm.stan`

file can be included immediately when initializing the package, any new Stan files can be added later to the `inst/stan`

folder. If we do not mind having `rstantools`

as a package dependency, it makes sense to set `auto_config = TRUE`

(the default), which avoids the need to manually reconfigure the package with `rstantools::rstan_config()`

whenever a `.stan`

file is `inst/stan`

are added, removed or modified.

## initialize R-package rstantools::rstan_create_package( path = "shinyStanModels", stan_files = "lm.stan" )

After updating the DESCRIPTION file and roxygenizing the package with `roxygen2::roxygenize()`

or `devtools::document()`

, the R-package can be installed with a call to `R CMD INSTALL`

or `devtools::install()`

. Note that building the package from source takes a while, since this is the moment when the Stan models are compiled and made available to R. The compiled Stan model originating from `lm.stan`

is now directly available in the internal object `stanmodels`

, a named list of S4-objects of class `"stanmodel"`

, with each S4-object containing the compiled model version of a single `.stan`

file in the `inst/stan`

folder:

shinyStanModels:::stanmodels[["lm"]] #> S4 class stanmodel 'lm' coded as follows: #> data { #> int<lower=0> N; #> vector[N] x; #> vector[N] y; #> } #> parameters { #> real alpha; #> real beta; #> real<lower=0> sigma; #> } #> model { #> y ~ normal(alpha + x * beta, sigma); #> } class(shinyStanModels:::stanmodels[["lm"]]) #> [1] "stanmodel" #> attr(,"package") #> [1] "rstan"

At this point we can directly sample from the S4-model objects with `rstan::sampling()`

:

system.time({ rstan::sampling( object = shinyStanModels:::stanmodels[["lm"]], data = list(N = 10L, x = seq_len(10), y = rnorm(10, seq_len(10), 0.1)), chains = 2, iter = 1000 ) }) #> #> SAMPLING FOR MODEL 'lm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 7e-06 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 0.024562 seconds (Warm-up) #> Chain 1: 0.013663 seconds (Sampling) #> Chain 1: 0.038225 seconds (Total) #> Chain 1: #> #> SAMPLING FOR MODEL 'lm' NOW (CHAIN 2). #> Chain 2: #> Chain 2: Gradient evaluation took 3e-06 seconds #> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds. #> Chain 2: Adjust your expectations accordingly! #> Chain 2: #> Chain 2: #> Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 2: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 2: #> Chain 2: Elapsed Time: 0.01885 seconds (Warm-up) #> Chain 2: 0.01167 seconds (Sampling) #> Chain 2: 0.03052 seconds (Total) #> Chain 2: #> user system elapsed #> 0.083 0.004 0.087

To keep everything together, we can just as well add the `app.R`

file to the R-package in a folder `inst/app`

. The contents of the `app.R`

file are now similar to our initial attempt, but with the call to `rstan::stan_model()`

left out:

library(shiny) library(rstan) ui <- fluidPage( sidebarLayout( sidebarPanel( numericInput("N", label = "N", value = 10) ), mainPanel( plotOutput("posteriors") ) ) ) server <- function(input, output, session) { ## draw samples draws <- reactive({ N <- input$N sampling( object = shinyStanModels:::stanmodels[["lm"]], data = list(N = N, x = seq_len(N), y = rnorm(N, seq_len(N), 0.1)), chains = 2, iter = 1000 ) }) ## plot histograms output$posteriors <- renderPlot({ req(draws()) op <- par(mfrow = c(1, 2), cex = 1.25) hist(extract(draws(), "alpha")[[1]], main = bquote("Posterior samples"~alpha), xlab = expression(alpha)) hist(extract(draws(), "beta")[[1]], main = bquote("Posterior samples"~beta), xlab = expression(beta)) par(op) }) } shinyApp(ui = ui, server = server)

The responsiveness of the shiny-app is now the same as in the previous section with the use the `rstanarm`

-package, but we are no longer constrained to only `rstanarm`

’s collection of Stan models.

## Models created with `brms`

Besides `rstanarm`

, the `brms`

-package also provides a flexible interface to build Stan models directly using R syntax. The difference between `rstanarm`

and `brms`

, however, is that `brms`

does not rely on pre-compiled Stan models and compiles generated `.stan`

files on-the-fly. This provides additional flexibility with respect to `rstanarm`

, but also means that calling `brms::brm()`

directly in an interactive shiny-app suffers from the same unresponsiveness as `rstan::stan_model()`

.

As a workaround, we can call `brms::make_stancode()`

to return the Stan program generated by `brms`

:

brms::make_stancode( formula = y ~ x, data = data.frame(x = numeric(1), y = numeric(1)), family = "gaussian" ) #> // generated with brms 2.14.4 #> functions { #> } #> data { #> int<lower=1> N; // total number of observations #> vector[N] Y; // response variable #> int<lower=1> K; // number of population-level effects #> matrix[N, K] X; // population-level design matrix #> int prior_only; // should the likelihood be ignored? #> } #> transformed data { #> int Kc = K - 1; #> matrix[N, Kc] Xc; // centered version of X without an intercept #> vector[Kc] means_X; // column means of X before centering #> for (i in 2:K) { #> means_X[i - 1] = mean(X[, i]); #> Xc[, i - 1] = X[, i] - means_X[i - 1]; #> } #> } #> parameters { #> vector[Kc] b; // population-level effects #> real Intercept; // temporary intercept for centered predictors #> real<lower=0> sigma; // residual SD #> } #> transformed parameters { #> } #> model { #> // likelihood including all constants #> if (!prior_only) { #> target += normal_id_glm_lpdf(Y | Xc, Intercept, b, sigma); #> } #> // priors including all constants #> target += student_t_lpdf(Intercept | 3, 0, 2.5); #> target += student_t_lpdf(sigma | 3, 0, 2.5) #> - 1 * student_t_lccdf(0 | 3, 0, 2.5); #> } #> generated quantities { #> // actual population-level intercept #> real b_Intercept = Intercept - dot_product(means_X, b); #> }

By including this Stan code in the `inst/stan`

folder and rebuilding the R-package, we circumvent the compilation step in `brms::brm()`

and can directly sample from the compiled Stan model with `rstan::sampling()`

as in the previous section. Note that the model input is slightly different, since `brms`

has generated the Stan code for a more general multiple linear model:

system.time({ brms_fit <- rstan::sampling( object = shinyStanModels:::stanmodels[["brms_lm"]], data = list(N = 10L, ## number of observations Y = rnorm(10, seq_len(10), 0.1), ## response vector K = 2L, ## number of predictors X = cbind(alpha = rep(1, 10), beta = seq_len(10)), ## predictor matrix prior_only = FALSE ## set to TRUE to evaluate only the priors ), chains = 2, iter = 1000 ) }) #> #> SAMPLING FOR MODEL 'brms_lm' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 6e-06 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 0.010256 seconds (Warm-up) #> Chain 1: 0.007216 seconds (Sampling) #> Chain 1: 0.017472 seconds (Total) #> Chain 1: #> #> SAMPLING FOR MODEL 'brms_lm' NOW (CHAIN 2). #> Chain 2: #> Chain 2: Gradient evaluation took 1e-05 seconds #> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds. #> Chain 2: Adjust your expectations accordingly! #> Chain 2: #> Chain 2: #> Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 2: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 2: #> Chain 2: Elapsed Time: 0.010079 seconds (Warm-up) #> Chain 2: 0.00762 seconds (Sampling) #> Chain 2: 0.017699 seconds (Total) #> Chain 2: #> user system elapsed #> 0.256 0.000 0.256 ## alpha is b_Intercept ## beta is b[1] rstan::summary(brms_fit, pars = c("b_Intercept", "b", "sigma"))[["summary"]] #> mean se_mean sd 2.5% 25% 50% #> b_Intercept 0.1362800 0.0025833278 0.07209531 -0.01553664 0.09773925 0.13826663 #> b[1] 0.9790349 0.0004261517 0.01127780 0.95785082 0.97213418 0.97865059 #> sigma 0.1014935 0.0019946898 0.03384615 0.05923768 0.07988739 0.09474943 #> 75% 97.5% n_eff Rhat #> b_Intercept 0.1801695 0.2614545 778.8522 1.000521 #> b[1] 0.9853540 1.0031303 700.3586 1.001592 #> sigma 0.1133669 0.1869625 287.9173 1.005812

**Remark:** the Stan code generated by `brms`

contains a bit of unnecessary complexity to sample from a simple linear model, but it does not take any effort to generate this Stan code, as we only need to provide the correct `brms`

model syntax.

## A note on deployment

When deploying the shiny-app to e.g. shinyapps.io or RStudio Connect using the `rsconnect`

-package, the R-package generated with `rstantools`

can be made available in a git repository on e.g. github or some other public repository, from which `rsconnect`

is able to fetch and install the R-package when deploying the shiny-app.

Another solution to deploy and host shiny-apps on a server is ShinyProxy, which launches shiny-apps from individual Docker containers. By installing the R-package generated by `rstantools`

when building the Docker image of the shiny-app, we ensure that we can directly sample from our compiled Stan models whenever a new Docker container is started. The following Dockerfile provides a minimal template to install the `rstantools`

-generated R-package from a bundled (`.tar.gz`

) package file and serve the shiny-app at port 3838:

FROM rocker/r-ver:4.0.3 # install system dependencies RUN apt-get update && \ apt-get install -y --no-install-recommends \ libv8-dev && \ apt-get clean && \ rm -rf /var/lib/apt/lists/ # install R packages (using littler) # this assumes .tar.gz exists in same folder as Dockerfile COPY shinyStanModels_0.1.tar.gz ./ RUN install2.r --error shiny rstan rstantools && \ install2.r --error shinyStanModels_0.1.tar.gz && \ rm *.tar.gz EXPOSE 3838 CMD ["R", "-e", "shiny::runApp(appDir = system.file('app', package = 'shinyStanModels'), port = 3838, host = '0.0.0.0')"]

For simplicity

`stan_glm()`

is used instead of`stan_lm()`

, as`stan_glm()`

automatically assigns*weakly informative priors*, whereas`stan_lm()`

expects a`prior`

argument using an additional call to`R2()`

.↩︎

**leave a comment**for the author, please follow the link and comment on their blog:

**R-bloggers | A Random Walk**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.